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FOREWORD 

Investigations  of  polynomial  operations  have  been  made  in  this  laboratory  for  a 
number  of  years.  Subroutines  have  been  prepared  on  the  basis  of  several  methods  of 
computation.  It  is  the  purpose  of  this  report  to  provide  analysis  and  documentation 
of  the  best  subroutines.  The  manuscript  was  completed  by  27  October  1977. 
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ABSTRACT 

Power  polynomials  and  trigonometric  polynomials  are  used  for  the  approximation 
of  general  functions.  Analysis  and  documentation  are  given  for  a  set  of  subroutines 
which  perform  polynomial  arithmetic,  find  roots  of  a  function,  or  construct  contours 
of  constant  function. 
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INTRODUCTION 

There  was  a  library  of  subroutines  on  the  Naval  Ordnance  Research  Calculator.  The 
subroutines  were  written  in  machine  language.  The  library  no  longer  was  useful  after 
NORC  was  dismantled.  Some  of  the  subroutines  have  been  transcribed  into  FORTRAN, 
and  are  described  herewith.  In  the  meantime,  the  NORC  library  has  been  eclipsed  by 
giant  library  projects.  There  are  the  System/360  Scientific  Subroutine  Package  of  the 
International  Business  Machines  Corporation1,  and  the  Subroutine  Library  of  the 
International  Mathematical  and  Statistical  Libraries,  Inc.8  Nevertheless,  the  subroutines 
from  the  NORC  library  present  capabilities  which  are  not  exactly  covered  by  the  giant 
libraries. 

Power  polynomials  and  trigonometric  polynomials  are  used  extensively  in  numerical 
analysis  because  they  can  be  evaluated  with  efficiency  on  digital  computers.  They  are 
used  often  to  approximate  other  functions  which  are  more  expensive  to  compute. 
They  give  continuous  interpolations  of  a  function  which  is  specified  only  at  discrete 
values  of  its  argument.  The  function  is  illustrated  graphically  as  a  set  of  ordinates  at 
a  set  of  abscissae.  The  function  can  be  represented  as  a  set  of  coefficients  of  the 
powers  of  its  argument  or  of  the  trigonometric  functions  of  its  argument.  Power 
polynomial  expansions  are  available  through  Lagrange  polynomial  approximation  or 
orthonormal  polynomial  approximation  for  arbitrary  arrangements  of  the  abscissae. 
The  power  polynomial  expansions  can  be  differentiated,  evaluated,  or  integrated. 

Lagrange  polynomials  are  convenient  because  the  polynomial  approximation  of  an 
arbitrary  function  is  given  directly  by  the  products  of  the  polynomials  and  discrete 
values  of  the  function.  In  the  construction  of  the  polynomials,  binomials  are  multiplied 
in  succession  to  form  a  polynomial  which  has  all  the  abscissae  for  roots,  then  this 
polynomial  is  divided  in  succession  by  binomials  to  obtain  the  Lagrange  polynomials. 
The  binomials  may  be  symbolic  insofar  as  they  are  represented  by  linear  functions 
of  their  argument.  There  is  a  loss  of  accuracy  to  rounding  error  in  the  derivation  of 
polynomials  of  high  order.  The  binomials  may  be  numeric  insofar  as  they  are  represented 
by  their  values  for  a  specific  value  of  their  argument.  There  is  no  significant  loss  of 
accuracy  to  rounding  error  in  the  evaluation  of  polynomials  of  any  order. 

The  orthonormal  polynomials  are  convenient  because  they  provide  a  least  squares 
approximation  for  a  function  when  the  number  of  discrete  values  of  the  function  is 
greater  than  the  order  of  the  approximation.  The  polynomials  may  be  orthonormal 
with  respect  to  integration  between  two  limits  of  integration,  or  they  may  be  orthonormal 
with  respect  to  summation  over  a  discrete  set  of  arguments.  In  the  construction  of 
the  polynomials  which  are  orthonormal  with  respect  to  summation,  a  matrix  of  values 
of  the  polynomials  is  synthesized  by  orthogonalization  and  the  coefficients  of  a 
three-term  recurrence  are  derived  by  summation.  The  coefficients  of  recurrence  may 
be  used  in  the  expansion  of  the  polynomials  into  representations  as  coefficients  of 
the  powers  of  the  argument,  or  the  coefficients  of  recurrence  may  be  used  in  the 
evaluation  of  the  polynomials  for  specific  values  of  the  argument. 

Many  high  -order  polynomials  have  such  large  coefficients  that  they  can  be  computed 
accurately  only  at  small  values  of  their  arguments.  Otherwise  the  rounding  error  in 
the  largest  term  of  their  expansion  exceeds  the  values  of  the  polynomials.  Expansion 
of  polynomials  into  representations  by  coefficients  is  feasible  only  for  polynomials  of 
low  degree. 
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A  high -order  polynomial  may  be  represented  as  a  set  of  discrete  values  with  an 
algorithm  for  interpolation,  or  as  a  recurrence  formulation  with  a  set  of  recurrence 
coefficients.  The  interpolation  and  the  recurrence  can  be  differentiated  or  evaluated 
but  not  integrated. 

High -accuracy  algorithms  for  integration  to  any  order  are  available  for  special 
spacing  between  abscissae.  Of  especial  importance  are  equal  spacing,  Chebyshev  spacing, 
and  Gauss  spacing. 

Equal  spacing  has  the  disadvantage  that  Lagrange  polynomials  of  high  degree  have 
large  amplitudes  of  oscillation  near  the  ends  of  a  set  of  abscissae.  Small  errors  in 
the  values  of  a  function  near  the  center  of  the  set  of  abscissae  cause  violent  oscillations 
near  the  ends. 

In  Chebyshev  spacing  the  abscissae  are  proportional  to  the  cosines  of  equally  spaced 
angles.  The  orthonormal  polynomials  for  Chebyshev  spacing  are  nearly  equal  to  the 
Chebyshev  polynomials.  The  amplitude  of  oscillation  of  the  polynomials  is  nearly 
uniform. 

In  Gauss  spacing  the  abscissae  are  proportional  to  the  roots  of  a  Legendre  polynomial. 
The  orthonormal  polynomials  for  Gauss  spacing  are  nearly  equal  to  the  Legendre 
polynomials.  The  order  of  accuracy  of  the  integration  of  a  function  is  nearly  twice 
the  number  of  values  of  the  function. 

In  the  discrete  Fourier  transform  the  orthogonality  of  trigonometric  functions  with 
respect  to  summation  is  used  to  obtain  the  coefficients  of  a  Fourier  series.  If  the 
order  of  the  Fourier  series  is  a  power  of  2,  then  there  is  available  the  fast  Fourier 
transform8.  Two  versions  of  the  FFT  are  known  as  the  Cooley-Tukey7  algorithm  and 
the  Sande-Tukey*  algorithm. 

Various  tactics  for  finding  roots  have  been  proposed  in  the  literature.  It  is  beyond 
the  scope  of  the  present  investigation  to  provide  a  literature  survey.  Libraries  of 
subroutines  have  been  assembled  by  various  organizations.  The  presence  of  an  algorithm 
in  such  libraries  is  evidence  that  the  algorithm  is  worthy  of  consideration. 

If  a  function  is  irregular  it  may  have  roots  anywhere.  A  set  of  values  of  the  function 
must  be  computed  at  a  closely  spaced  set  of  abscissae.  The  values  are  sorted  in  order 
to  determine  where  the  function  changes  sign. 

If  a  function  is  represented  by  a  smooth  curve  of  ordinates,  then  the  presence  of 
roots  can  be  predicted  from  the  behavior  of  the  function  between  roots.  Only  a  few 
computations  of  ordinates,  slopes,  and  curvatures  are  required  in  order  to  locate  and 
determine  the  roots. 

In  the  method  of  false  position  or  regula  falsi3'4,  two  abscissae  which  straddle  a 
root  must  first  be  found.  The  values  of  the  function  are  opposite  in  sign  at  the  two 
abscissae.  Linear  interpolation  is  used  to  obtain  a  new  abscissa.  The  interpolation  is 
repeated  between  the  nearest  abscissae  for  which  the  function  has  opposite  signs.  The 
curvature  of  the  curve  of  ordinates  tends  to  make  the  sign  of  the  function  at  each 
new  abscissa  uniformly  the  same  and  opposite  to  the  sign  of  the  function  at  a  fixed 
abscissa.  Although  the  abscissae  improve  gradually  in  accuracy,  the  rate  of  convergence 
is  geometric,  because  each  new  abscissa  is  obtained  with  a  slope  which  does  not 
continue  to  improve  in  accuracy. 

In  the  Pegasus18  method,  the  location  of  the  new  abscissa  overshoots  the  regula 
falsi  by  an  amount  which  gives  improved  convergence. 

In  the  Lehmer  Schur19  method,  the  field  is  partitioned  into  areas  by  an  arrangement 
of  circles  with  various  centers  and  radii.  The  coefficients  of  polynomials  with  complex 
argument  are  combined  with  the  coefficients  of  polynomials  with  inverse  argument  so 
as  to  eliminate  the  terms  of  highest  degree.  When  the  polynomials  have  been  reduced 
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to  a  constant,  then  the  sign  of  the  constant  indicates  whether  a  root  is  located  within 
the  circle  of  inversion.  Repetition  with  the  various  circles  identifies  the  areas  within 
which  the  roots  are  located.  The  disadvantage  of  the  method  is  that  the  time  for 
reduction  of  polynomials  for  each  circle  is  proportional  to  the  square  of  the  order, 
and  the  method  is  too  expensive. 

In  the  Newton- Raphson3,4  method,  a  value  of  the  function  and  its  derivative  at  an 
iterant  is  used.  At  the  iterant  a  tangent  to  the  curve  of  ordinates  is  constructed.  The 
point  of  intersection  of  the  tangent  with  the  axis  of  abscissae  is  the  new  iterant. 

In  the  secant  method,  the  values  of  the  function  at  two  successive  iterants  are 
used.  A  chord  is  constructed  between  the  ordinates  at  the  two  iterants.  The  point  of 
intersection  of  the  chord  with  the  axis  of  abscissae  is  the  new  iterant. 

The  advantage  of  the  tangent  is  that  it  does  not  require  as  many  iterations.  The 
advantage  of  the  chord  is  that  it  does  not  require  the  differentiation  of  the  function. 
The  rates  of  convergence  are  quadratic,  because  the  accuracy  of  the  abscissae  and 
the  accuracy  of  the  slope  both  improve  together. 

In  the  Muller  method20,  the  values  of  the  function  at  three  successive  iterants  are 
used.  The  three  values  of  the  function  are  fitted  with  a  quadratic  polynomial.  The 
nearer  root  of  the  quadratic  polynomial  is  the  next  iterant.  Convergence  is  complete 
if  the  increment  from  the  iterant  to  the  Newton-  Raphson  intercept  is  less  than  a 
tolerance.  The  three  consecutive  iterants  which  define  the  quadratic  polynomial  must 
be  distinct.  Otherwise  the  roots  of  the  quadratic  polynomial  may  be  out  of  bounds. 
In  that  case,  the  iterant  is  brought  in  to  a  circle  with  radius  equal  to  the  geometric 
mean  of  the  distances  to  adjacent  roots. 

In  the  method  of  Jarrett  and  Nudds21,  the  quadratic  polynomial  is  replaced  by  a 
linear  rational  function. 

In  the  method  of  Jenkins  and  Traub22,  the  function  is  a  polynomial.  A  set  of  auxiliary 
polynomials  is  constructed  by  a  recurrence  relation.  A  fraction  of  the  function  is 
subtracted  from  each  auxiliary  polynomial  to  give  a  difference  which  is  zero  at  an 
iterant.  The  difference  is  divided  by  a  binomial  to  give  auxiliary  polynomials  of  one 
lower  degree  than  the  function.  At  all  roots  of  the  function  except  the  nearest  root, 
the  values  of  the  auxiliary  polynomials  approach  zero  by  a  powering  process. 
Newton -Raphson  iteration  is  applied  to  the  quotients  between  the  function  and  the 
auxiliary  polynomials.  Insofar  as  all  roots  of  the  auxiliary  polynomials  except  one 
approach  the  roots  of  the  function,  the  Newton-Raphson  iteration  is  applied  to  a 
progressively  more  nearly  linear  function. 

In  the  present  system,  it  is  assumed  that  the  best  tactic  is  to  go  directly  in  that 
direction  in  which  a  root  is  most  likely  to  lie.  There  are  two  stages  in  the  search  for 
roots.  In  the  hunting  stage,  the  function  is  approximated  locally  at  an  iterant  by  a 
quadratic  polynomial.  Steps  of  predetermined  length  are  taken  away  from  the  iterant 
in  the  direction  of  the  nearer  root  of  the  quadratic  polynomial.  Convergence  of  the 
Newton  Raphson  iteration  is  tested  with  the  aid  of  a  criterion.  The  hunting  stage 
terminates  when  a  region  of  convergence  has  been  reached.  In  the  homing  stage,  the 
intercept  of  the  tangent  with  the  complex  plane  is  computed  at  each  iterant.  The 
intercept  becomes  the  new  iterant  as  long  as  the  displacement  of  the  intercept  is  less 
than  the  displacement  of  the  iterant.  Termination  of  the  homing  stage  occurs  if  the 
displacement  of  the  intercept  is  equal  to  or  greater  than  the  displacement  of  the 
iterant.  Termination  may  occur  because  rounding  error  has  upset  the  convergence  or 
because  the  iterant  has  moved  out  of  the  region  of  convergence.  Whether  termination 
is  caused  by  rounding  error  or  by  nonconvergence  is  tested  by  comparison  between 
the  displacements  and  a  tolerance.  If  the  displacements  are  less  than  the  tolerance, 
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the  iterant  is  accepted  as  a  root.  If  the  displacements  are  more  than  the  tolerance, 
the  hunting  stage  is  resumed  with  double  the  tolerance. 

After  each  root  of  a  polynomial  has  been  determined,  the  root  of  the  polynomial 
may  be  eliminated  by  a  synthetic  division  to  obtain  a  reduced  polynomial  of  one  lower 
degree.  Experience  has  shown  that  with  each  reduction  of  the  degree  of  the  polynomial 
the  remaining  roots  are  disturbed  and  the  accuracy  deteriorates.  Reduction  of  degree 
is  not  feasible  for  transcend  dal  functions.  In  the  present  system  the  function  is 
replaced  by  the  quotient  between  the  function  and  that  polynomial  whose  roots  are 
equal  to  all  roots  already  determined. 

Contours  are  used  for  the  illustration  of  information  which  may  range  from  irregular 
data  to  smooth  functions.  The  contours  for  irregular  data  are  based  upon  local 
interpolations  within  one  element  of  a  grid,  whereas  contours  for  smooth  functions 
can  be  based  upon  global  formulations  for  the  entire  field. 

An  equation  between  two  coordinates  defines  implicitly  a  relationship  which  is 
represented  in  a  graph  of  the  coordinates  by  a  contour  line.  The  equation  is  solved 
by  an  iteration  which  starts  with  trial  values  of  the  coordinates  and  follows  a  hunting 
procedure  and  a  homing  procedure.  During  hunting  the  trial  point  is  displaced  by 
steps  of  predetermined  length.  Before  the  contour  is  reached  the  steps  are  perpendicular 
to  the  contour  while  after  the  contour  is  reached  the  steps  are  parallel  to  the  contour. 
After  a  contour  is  reached  it  is  followed  until  it  goes  out  of  bounds  or  closes  on  itself. 
If  it  does  not  close  on  itself,  then  the  trial  point  returns  to  the  initial  point  and  the 
other  half  of  the  contour  is  followed  until  it  also  goes  out  of  bounds. 

Two  programs  on  NORC  were  designed  to  trace  section  lines  on  a  ship  or  contour 
lines  in  a  LORAN  map.  These  programs  were  prepared  when  the  available  cathode-ray 
printers  were  dot  plotters.  The  rasters  were  coarse  and  the  simulation  of  curves 
required  care  in  the  arrangement  of  dots.  Later  the  rasters  were  finer  and  the 
cathode  ray  printers  were  vector  plotters.  Now  the  rasters  are  still  finer  and  the 
latest  cathode -ray  printers  again  are  dot  plotters.  The  dots  are  arranged  so  as  to 
simulate  vectors.  Any  simulation  of  a  smooth  curve  by  a  vector  plotter  necessarily  is 
a  polygonization.  The  simulation  of  a  smooth  curve  by  a  polygonization  requires  an 
optimization  of  the  lengths  of  the  sides  of  the  polygonization. 


POLYNOMIAL  ARITHMETIC 

Analysis 

Let  the  polynomials  A(z),  B(z),  C(z)  be  defined  by  the  equations 


L-l 

A(z)  E  akz* 

A/-1 

(i) 

B(z)  -  E  <>kzk 

fc= 0 

JV-l 

(2) 

r(z)  =  E  ckzk 

k-0 

(3) 

where  ak,bk.ck  are  real  coefficients  and  L.M.N  are  the  numbers  of  coefficients 
respectively  The  argument  z  may  be  real  or  complex.  The  coefficients  of  each  polynomial 
are  stored  in  an  array  in  ascending  order.  If  the  degree  of  the  polynomial  is  less  than 
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the  size  of  its  array,  then  the  upper  end  of  the  array  is  filled  with  zeros.  The  dimensions 
of  an  array  are  given  by  a  two-integer  specification.  The  first  integer  is  the  interval 
between  elements  of  the  array,  and  the  second  integer  is  the  number  of  elements  in 
the  array. 

In  the  transfer  of  a  polynomial,  the  coefficients  of  the  new  polynomial  are  related 
to  the  coefficients  of  the  old  polynomial  in  accordance  with  the  equation 

i>*  =  a*  (4) 

Transfer  continues  until  the  second  array  has  been  filled  with  coefficients  from  the 
first  array,  or  with  zeros  if  the  first  array  is  shorter  than  the  second  array. 

In  the  addition  of  two  polynomials,  the  coefficients  of  the  new  polynomial  are  related 
to  the  coefficients  of  the  old  polynomials  in  accordance  with  the  equation 

c*  =  a*  +  bk  (5) 

Addition  continues  until  the  new  array  is  filled. 

In  the  subtraction  of  two  polynomials,  the  coefficients  of  the  new  polynomial  are 
related  to  the  coefficients  of  the  old  polynomials  in  accordance  with  the  equation 

c*  =  -  bk  (6) 

Subtraction  continues  until  the  new  array  is  filled. 

In  the  multiplication  of  two  polynomials,  the  coefficients  of  the  new  polynomial  are 
related  to  the  coefficients  of  the  old  polynomials  in  accordance  with  the  equation 

Cm  =  S  O-m-kblc  ( 7 ) 

k=0 

The  indices  to  k  and  k  are  coordinates  in  a  rectangle  with  dimensions  equal  to  the 
lengths  of  the  old  polynomials.  In  the  summation  the  indices  advance  diagonally  across 
the  rectangle.  Initialization  of  the  summation  is  at  the  left  or  bottom  of  the  rectangle, 
and  termination  of  the  summation  is  at  the  top  or  right  of  the  rectangle. 

In  the  division  of  polynomials,  fractions  of  the  denominator  are  subtracted  from 
the  numerator  until  the  numerator  has  been  eliminated.  The  fractions  which  are 
applied  to  the  denominator  are  the  coefficients  in  the  quotient.  Elimination  is 
accomplished  with  the  aid  of  the  substitutions 


^0 


(8) 


am*k  "*  am+k  cm^k  (®) 

where  k  is  limited  to  the  range  k  <  M  and  to  the  range  k  +  m  <  ,V. 

The  quotient  between  two  polynomials  is  a  finite  polynomial  if  all  roots  of  the 
denominator  coincide  with  roots  of  the  numerator.  Otherwise  the  quotient  is  an  infinite 
series.  The  circle  of  convergence  of  the  infinite  series  goes  through  that  nearest  root 
of  the  denominator  which  does  not  correspond  to  a  root  of  the  numerator. 

Polynomial  arithmetic  with  real  coefficients  can  be  adapted  to  polynomial  arithmetic 
with  complex  coefficients.  Let  complex  polynomials  be  given  by  the  expressions 

A  +  iB  C+iD  (10) 


Then  complex  transfer  is  expressed  by  the  substitutions 

C  -*A  D -B 


(11) 
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Complex  addition  is  expressed  by  the  equation 

(A  +  iB)  +  (C  -f-  iD)  =  (.4  +  C)  +  ( B  +  B)i 
Complex  subtraction  is  expressed  by  the  equation 

(A+iB)-(C  +  iD)  =  (A-C)  +  (B-D)i  (13) 

Complex  multiplication  is  expressed  by  the  equation 

(A  +  iB)(C  +  iD)  =  (AC  -  BD)  +  (AD  +  BC)i  (14) 

Complex  division  is  expressed  by  the  equation 

A  +  iB  (AC  +  BD)  +  (BC  -  AD)i 

C+iD~  CC  +  DD  (15^ 

The  polynomials  A,  B.  C,  D  have  real  coefficients  but  complex  arguments. 
Programming 

SUBROUTINE  PLNM  (MO.  PA,  NA,  PB,  NB,  PC,  NC) 

*************************+************************+****************************************** 
FORTRAN  SUBROUTINE  FOR  POLYNOMIAL  ARITHMETIC 


The  mode  of  operation  is  given  in  argument  MO.  The  polynomials  A,B,  C  are  given  or 
stored  in  the  arrays  PA,  PB,  PC.  The  specifications  of  the  polynomials  are  given  in  the 
arrays  NA,  NB,  NC.  The  specification  for  a  polynomial  PX  consists  of  the  array  NX  as 
given  in  the  following  table: 

Address  Specification 

NX(1)  Interval  between  coefficients 

NX(2)  Number  of  coefficients 


There  are  no  compatibility  requirements  on  the  lengths  of  the  polynomials.  The 
repertory  of  operations  and  call  lines  is  given  in  the  following  table: 


Operation 

Call 

Line 

B  -  A 

CALL 

PLNM 

C  A  +  B 

CALL 

PLNM 

C  A  -  B 

CALL 

PLNM 

C  AB 

CALL 

PLNM 

CAB 

CALL 

PLNM 

(1,  PA,  NA,  PB,  NB) 

(2,  PA,  NA,  PB,  NB,  PC,  NC) 
(3,  PA,  NA,  PB,  NB,  PC,  NC) 
(4,  PA,  NA,  PB,  NB,  PC,  NC) 
(5.  PA.  NA,  PB,  NB,  PC,  NC) 


During  polynomial  division,  leading  zeros  are  shifted  off  to  the  left  as  far  as  possible 
in  the  numerator  and  the  denominator. 


MULTIPLEX  POLYNOMIAL  EVALUATION 


Analysis 

Let  ip(u)  be  a  polynomial  of  (n  -  l)th  degree  in  the  argument  u.  Let  the  polynomial 


be  expressed  in  terms  of  its  argument  by  the  equation 

u- 1 

<p(u)  =  £  akuk  (16) 

t=o 

where  the  coefficient  at  is  the  kth  element  of  an  array.  The  first  derivative  of  the 
polynomial  is  given  by  the  equation 

d<p(u)  n~l 

— t~  =  E  kakuk  (17) 

du  k=0 

and  the  second  derivative  of  the  polynomial  is  given  by  the  equation 

=  V  k{k  -  l)0tu*-*  (18) 

Jfc  —  0 

Integration  of  the  polynomial  is  expressed  by  the  equation 

f  <fi(u)du=  £  ut+‘  (19) 

Jo  k= 0  ^  +  1 

Whether  the  operation  is  integration,  evaluation,  first  differentiation,  or  second 
differentiation  is  determined  by  a  mode-  of-  operation  parameter. 

Programming 

SUBROUTINE  MPLNMV  (WO,  AU,  NC,  AC,  FV) 

*****tf^********************4^*********4^*******************4*****S*»************************* 

FORTRAN  SUBROUTINE  FOR  MULTIPLEX  POLYNOMIAL  EVALUATION 

***************-************.**.*******.***i***.***************.**.*m******m**’ti*-****».****.*^t**:***.**** 

The  mode  of  operation  is  integration  if  MO  =  - 1.  evaluation  if  WO  =  0,  first  differentiation 
if  MO  -  -*-l,  and  second  differentiation  if  WO  -  -^2.  The  argument  u  is  given  in  the 
address  AU.  The  number  of  coefficients  n  is  given  in  the  address  NC.  The  coefficients 
a0,  a„.  [  are  given  in  the  n- array  AC.  The  value  of  the  function  is  stored  in  address  FV. 


LAGRANGE  POLYNOMIAL  EXPANSION 

Analysis 

Let  ut,  ••',  un  be  discrete  values  of  the  argument  u,  and  let  ^,(w),  <pn{u)  be  Lagrange 
polynomials.  The  mth  polynomial  is  defined  by  the  equation 


.  t.n  _  ( u  -  ui)'(u  -  um  ,)(u  -  um+l)  -(u  -  un) 

/  x  /  w  x  I  X  (20) 

(\  -  «l)-K  -  Wm_,)(ltm  -  -  Un) 

The  Lagrange  polynomials  have  the  property  which  is  expressed  by  the  equation 

=  (21) 

where  6mk  is  zero  if  k  r.  m  but  is  unity  if  k  =  m.  The  Lagrange  polynomials  are  expressed 
in  terms  of  powers  of  the  argument  by  the  equation 

TV-  I 

<Pm(u)  =  E  V1  (22) 
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T 


where  the  coefficient  is  the  mfeth  element  of  a  matrix  of  coefficients. 

Let  the  partial  polynomial  Fm(u)  be  defined  by  the  equation 

Fm  (u)  =  {u-ul)-(u-um)  (23) 

and  be  expressed  in  terms  of  powers  of  the  argument  by  the  equation 

m 

f'mM  =  E  (24) 

fc=0 

The  coefficients  of  the  partial  polynomials  are  generated  with  the  aid  of  the  recurrence 
equations 


Cmk  C m-l.k-l  (25) 

and 

^mm  —  1  (26) 

The  coefficients  of  the  partial  polynomials  are  stored  temporarily  in  the  matrix  of 
coefficients. 

Cycling  of  the  recurrence  leads  to  the  complete  polynomial  F(u)  which  is  defined 
by  the  equation 


F(u)  =  (u- ut)  -(u-un)  (27) 

and  is  expressed  in  terms  of  powers  of  the  argument  by  the  equation 

F(u)  =  £  ctuk  (2B) 

*=o 

The  coefficients  of  the  complete  polynomial  are  stored  in  the  (n  +  l)th  row  of  the 
matrix  of  coefficients. 

The  Lagrange  polynomial  ^m(it)  is  recovered  from  the  complete  polynomial  F(u) 
through  division  by  the  binomial  u  -  um.  The  algorithm  for  division  is  just  the  nested 
method  for  polynomial  evaluation.  The  coefficients  of  the  Lagrange  polynomials  are 
given  initially  by  the  recurrence  equation 

®m*  =  Cfi  +  (29) 

and  then  they  are  normalized  through  division  by  the  product 

(^•m  ^1 )  '(^m  —  ^m+1  )' "  (^An  ~  ^n)  (30) 

All  coefficients  are  stored  in  the  computer  in  ascending  order. 

Inasmuch  as  the  mth  Lagrange  polynomial  is  unity  at  the  mth  argument  and  zero 
for  other  discrete  arguments,  an  arbitrary  function  f(u)  is  expressed  by  the  equation 

n 

f{u)=  E  (31) 

m-~  1 

The  polynomial  approximation  of  an  arbitrary  function  is  derived  from  a  set  of  values 
of  the  function  at  the  discrete  arguments  through  vector-matrix  multiplication. 
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Programming 

SUBROUTINE  LGRNGX  (AU,  NA,  AC) 

ft******************************************************************************************** 

FORTRAN  SUBROUTINE  FOR  LAGRANGIAN  POLYNOMIAL  EXPANSION 
*******%**#*******************#**********************•*****************  <.********************* 

The  arguments  itj.---.Un  are  given  in  the  n- array  AU.  The  order  n  is  given  in  the 
argument  NA.  The  coefficients  of  the  polynomials  <fiy(u),  — .  <fn(u)  are  stored  in  the 
(n  +  1)  x  n  array  AC.  The  coefficients  of  the  polynomial  F{u)  are  stored  in  the  (n  +  l)th 
row  of  the  array  AC. 


LAGRANGE  POLYNOMIAL  EVALUATION 


Analysis 

If  u*uk  for  any  k,  then  the  mth  Lagrange  polynomial  y>m(u)  >s  evaluated  by  the 
equation 

,  ^  =  (u  -  Ui)-(u  -  um-t)(u  -  um+i)-ju  -  Un) 

("^■m  ~  )’" ("M-m  ~  ^m- 1 )("^m  ^m+l )‘ "(V*m  ~  ^-n) 

The  normalization  divisor  of  the  mth  polynomial  is  given  by  the  product 
(um  -  U,)  --(um  -  -  um+1)-  (um  -  -Un) 

which  is  computed  initially  from  the  discrete  data.  The  function  F(u)  in  the  equation 

/’(u)  =  (u-u1) --(u-u „)  (34) 


(32) 

(33) 


is  computed  from  the  given  argument  u,  is  divided  by  each  difference  u  -  um  and  by 
the  normalization  divisor  to  give  the  values  of  the  successive  polynomials  ^m(u). 

The  first  derivative  of  the  polynomial  <pm(u)  is  given  by  the  equation 


d^-{u)  =  + 


1 


1 


1 


u  ~  u. 


u  -  Um-1 


•u-iWi 


U-Un 


<pm(u)  (35) 


The  sums  of  the  reciprocals  of  the  first  m  -  1  differences  u  -  u^.  j  are  accumulated 
and  are  stored  in  an  array. 

The  second  derivative  of  the  polynomial  <pm{u)  is  given  by  the  equation 


du‘ 


1 


1 


1 


1 


(U-Uj)2 

1 


(u  -  Um_i)2 
1 


(u 


V-T, 

1 


u  -  «m-l 


(u  -  Un)Z 

1 _ 

u-un 


9>mM 


9>mM  06) 


Polynomial  expansion  and  cancellation  of  squares  of  reciprocals  of  differences  reduce 
the  second  derivative  to  the  summation  in  the  equation 


dz  n  n  1 

-r~z  <pm(u)  =  E  E  - w- - r 

du  i=u=i  (u-uJiu  -Uf) 


¥>m(u) 


(i*j)  (37) 


where  the  summation  does  not  include  any  terms  for  which  i  =  m,  j  =  m,  or  i  j. 

The  terms  in  the  summation  may  be  arranged  in  the  form  of  a  matrix  such  that 
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the  i.jth  term  occupies  the  i.jth  position  in  the  matrix.  The  mth  row,  the  mth  column, 
and  the  diagonal  of  the  matrix  contain  zeros.  There  is  a  matrix  for  each  value  of  m. 
The  sum  of  elements  in  the  upper  left  quadrant  is  twice  the  sum  of  products  of  the 
reciprocal  of  the  difference  u  -  and  the  sum  of  the  first  m  -  2  reciprocals  of  the 
differences  u  -  itm-a-  The  sums  of  the  elements  in  the  upper  left  quadrant  are  computed 
and  are  stored  in  an  array. 

The  computation  of  the  second  derivative  is  completed  with  further  accumulations 
of  the  products  of  the  reciprocals  of  differences.  The  sum  of  elements  in  the  lower 
right  quadrant  is  twice  the  sum  of  products  of  the  reciprocal  of  the  difference  u  -  um+ , 
and  the  sum  of  the  last  n  -  m  -  1  reciprocals  of  the  differences  u  -  um+1.  The  sum  of 
elements  in  the  lower  left  quadrant  and  in  the  upper  right  quadrant  is  twice  the 
product  of  the  sum  of  the  first  m  -  1  reciprocals  of  differences  and  the  sum  of  the 
last  n  -  m  reciprocals  of  differences. 

The  computation  of  the  first  derivative  is  completed  with  further  accumulations  of 
the  reciprocals  of  differences.  The  sum  of  the  first  m  -  1  reciprocals  of  differences  is 
added  to  the  sum  of  the  last  n  -  m  reciprocals  of  differences  to  complete  the  summation. 

If  u  =  uk  for  any  k,  then  the  reciprocal  of  the  difference  u  -  uk  is  omitted  throughout 
the  computation,  and  the  values  of  the  polynomials  are  replaced  by  zero  when 

m  *  k  and  by  unity  when  m  =  k. 

Programming 

SUBROUTINE  LGRNGN  (AU,  NA,  AN) 

************************************************+*****************♦*******************•****** 
FORTRAN  SUBROUTINE  FOR  LAGRANG1AN  NORMALIZATION  DIVISORS 
********************************************************************************************* 

The  arguments  u,,  ■,un  are  given  in  the  n- array  AU.  The  number  of  arguments  is 
given  in  the  argument  NA.  The  normalization  divisors  are  stored  in  the  n  array  AN. 

SUBROUTINE  IGRNGV  (MO,  NA,  QU,  AU,  AN,  FF,  DF.  SF) 

*************************************************«**********************«***•**+«*«***«•***** 
FORTRAN  SUBROUTINE  FOR  LAGRANGIAN  POLYNOMIAL  EVALUATION 

*************4r***********44**4*****4***4*4******************************»*»*«*«*********«***+ 

The  mode  of  operation  MO  is  0  for  functions,  1  for  first  derivatives,  and  2  for  second 
derivatives.  The  number  of  discrete  arguments  n  is  given  in  argument  NA.  The  variable 
argument  u  is  given  in  argument  OU.  The  discrete  arguments  u,,  ,un  arc  given  in’ 
the  n  array  AU,  and  the  normalization  divisors  are  given  in  the  n  array  AN.  The 
functions  are  stored  in  the  n  array  FF,  the  first  derivatives  are  stored  in  the  n  array 
DF,  and  the  second  derivatives  are  stored  in  the  n  array  Sf . 


ORTHONORMAL  POLYNOMIAL  SYNTHESIS 

Analysis 

Let  it,,  -  ,u„  be  the  discrete  values  of  the  argument  u  at  n  stations,  and  let 
<Po(u),  •••,  <pn  ,(u)  be  n  polynomials  of  progressively  increasing  degree.  The  polynomials 


are  orthonormal  if  they  satisfy  the  condition 

n 

£  (PiKlOik)  =  <5*m  (38) 

i=l 

where  5^  is  zero  if  fc  *  m  and  is  unity  if  k=m.  Any  polynomial  of  degree  m  can  be 
expanded  in  terms  of  the  orthonormal  polynomials  of  degrees  0  through  m.  The 
coefficients  in  the  expansion  are  given  by  the  solution  of  a  sequence  of  equations 
which  start  at  the  highest  degree  and  compare  terms  of  the  same  degree.  Let  <pm(u) 
be  expressed  by  the  equation 

TO-i 

=7m.U<Pm-iM  +  £  7v<PvM  (39) 

v=0 

where  the  y„  are  constant  coefficients.  Multiplication  throughout  by  <pk(u)  and 
summation  isolates  yk  in  accordance  with  the  equation 

n  n 

£  VkMVmM  =  7m  £  +  yk  =  0  (40) 

<= 1  i= 1 

The  polynomial  u<pk(u)  may  be  expressed  by  the  equation 

*+i 

u<fik(u)  =  £  ev<e„(u)  (41) 

v~0 

where  the  e„  are  constant  coefficients.  All  terms  of  the  summation  are  orthogonal  to 
(Bm-iW  and  the  coefficient  yk  therefore  is  zero  when  kSm-3.  The  orthonormal 
polynomials  can  be  generated  by  a  three-term  recurrence  equation. 

The  polynomial  of  zero  degree  is  given  by  the  equation 


Vo  M  ~ 


1 

Vn 


(42) 


and  the  polynomial  of  first  degree  is  given  by  the  equation 

Vi{u)  =  —  (it  -  Po)VoM  (43) 

«o 

where  the  constants  a0  and  /S0  are  defined  by  the  equations 

I  n 

a0  =  ~7=  £  ««?»(««)  (44) 

vn  i=i 


and 


Po  —  £ 


(45) 


The  subsequent  polynomials  are  generated  by  the  three-term  recurrence  equation 


~  ^m-it")  -  Pm-lVm-lM  ~  <*m-zVm-i(u)\ 


(46) 
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where  the  constants  nm_!  and  iff*,-,  are  defined  by  the  equations 

n 

am-i  =  T,ui<fim.,(ui)<pm(ui)  (47) 

t=l 

and 

n 

0m-l  =  E  ifaJ'Pm-lfal)  (48) 

1=1 

The  coefficients  am-],  0m-i  are  stored  in  a  row  array  and  the  values  ^(Ut)  are  stored 
in  a  matrix  array. 

Let  a  function  f(u)  be  expanded  as  a  series  in  the  orthonormal  polynomials  as 
expressed  by  the  equation 


f(u)  =  Y.  cm<fim(u)  (49) 

m=0 

The  square  of  the  standard  deviation  a 2  is  given  by  the  equation 

c2="  -  E  E  c^u*)!  (50) 

n  1=1  (  k=0  t 

Differentiation  with  respect  to  a  coefficient  gives  the  equation 

-  zr~  =  \  E  j/(«i)  -  E  ct^t(ut)|  ^(u*)  (51) 

8cm  n  <=1  (  t=0  ) 

For  least  squares  the  coefficients  are  given  by  the  equation 

n 

Cm  =  E  /(Wl)^m(Ml)  (52) 

i=l 

The  square  of  the  standard  deviation  is  given  by  the  equation 

|E/2(ui)-  E‘cmzj  (53) 

n  1  1=1  m=0  > 

For  a  complete  expansion  in  terms  of  all  n  orthonormal  polynomials  the  function  is 
equal  to  its  series  representation  at  every  discrete  argument  and  the  standard  deviation 
is  zero. 

Programming 

SUBROUTINE  ORTHOS  (AU,  MA,  AA,  NA,  AR) 

******^******1,************************************** ***********************************  ****** 

FORTRAN  SUBROUTINE  FOR  ORTHONORMAL  POLYNOMIAL  SYNTHESIS 
* ******************************************************************************************** 

The  arguments  ult  ••■,Tzn  are  given  in  the  1  xn  array  AU.  The  number  of  polynomials 
m  is  given  in  argument  MA.  The  values  of  the  polynomials  <p0(u),  •••,  (fftn_I(u.)  are  stored 
in  the  m  x  n  array  AA.  The  number  of  arguments  n  is  given  in  argument  NA.  The  recurrence 
coefficients  a0,  (S0,  •••,  am_z,  /Sm_2  are  stored  in  the  2m.  -  2  array  AR. 
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ORTHONORMAL  POLYNOMIAL  EXPANSION 

Analysis 

Let  the  mth  orthonormal  polynomial  ^m(u)  be  expressed  by  the  equation 

n- 1 

*>m(**)  =  £ 

*-0 

The  polynomial  of  zero  degree  is  given  by  the  equation 

?o(u)  =  -t= 

Vn 

for  which  all  coefficients  are  zero  except 


The  polynomial  of  first  degree  is  given  by  the  equation 

<Pi(u)  =  —  (u  -  p0)<fioM 
ao 


for  which  all  coefficients  are  zero  except 

J0 

al0  -  " 


Thereafter  the  polynomials  of  progressively  increasing  degree  are  given  by  the  equation 
?m(u)  -  1  -  Pm-iVm-lM  -  am_z<pmJu)\  (59) 

I  \  ' 

for  which  the  coefficients  are  generated  with  the  aid  of  the  recurrence  equation 

“mi  =  ~  ”  Pm-iUm-l,*  ~  (60) 

am  1  \  ) 

The  constants  am.t,  Pm  i  are  given  in  a  row  array  and  the  coefficients  are  stored 
in  a  matrix  array. 

If  a  function  /(u)  is  represented  by  a  series  of  the  orthonormal  polynomials,  then 
the  function  f(u)  is  converted  into  a  series  of  the  powers  of  the  argument  u  by  a 
vector- matrix  multiplication. 

Programming 

SUBROUTINE  ORTHOV  (NA,  AR,  NC.  AC) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  ORTHONORMAL  POLYNOMIAL  EXPANSION 
***********  ****************************************  ****************************  ************** 

The  number  of  arguments  n  is  given  in  the  argument  NA.  The  recurrence  coefficients 
ag,p0,  ■■■,  am.z,pm  2  are  given  in  the  2 m  -  2  array  AR.  The  number  of  coefficients  m  is 
given  in  the  argument  NC.  The  coefficients  of  the  polynomials  ?o(u),  •••,  pm_t(u)  are 
stored  in  the  m  x  m  array  AC. 
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ORTHONORMAL  POLYNOMIAL  EVALUATION 

Analysis 

The  polynomial  of  zero  degree  is  given  by  the  equation 

1 


niv)  = 


Vn 


and  the  polynomial  of  first  degree  is  given  by  the  equation 

?,(u)  =  —  (u-  p0)VoM 
«o 


(61) 


(62) 


The  polynomials  of  progressively  higher  degree  are  generated  by  the  three-term 
recurrence  equation 


(63) 


<fim(u)  =  — l—  ju^m_t(u)  -  -  am-2^m-a(u)( 

«m-l  (  > 

where  constants  am_i  and  are  given  in  a  2m-  2  array. 

The  first  derivative  of  the  polynomial  of  zero  degree  is  given  by  the  equation 

<Pq(u)  =  0  (64) 

du 


and  the  first  derivative  of  the  polynomial  of  first  degree  is  given  by  the  equation 


ViM  =  ~  ?o(«) 

du  a0 


(65) 


The  first  derivatives  of  the  polynomials  of  progressively  higher  degree  are  generated 
by  the  three-term  recurrence  equation 


d 

du 


<Pm(u)  =  — —  Um-lM  +  U<Pm-l(u)  -  ~ 

am-l  \  ' 


(66) 


where  ^(u)  is  the  first  derivative  of  < pm(u). 

The  second  derivative  of  the  polynomial  of  zero  degree  is  given  by  the  equation 


du‘ 


<p0(u)  =  0 


(67) 


and  the  second  derivative  of  the  polynomial  of  first  degree  is  given  by  the  equation 

d' 


du' 


Vi  (“)  =  0 


(68) 


The  second  derivatives  of  the  polynomials  of  progressively  higher  degree  are  generated 
by  the  three-term  recurrence  equation 

~g  VmM  =  Upm-tM  +  Wm-lM  -  Pm-lPm-lM  ~  «m-g? m-g(u)j 
du“  «m  -l  (  > 

where  >s  the  second  derivative  of  <pm(u). 


(69) 
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Programming 

SUBROUTINE  ORTHOV  (MO.  NA.  AU.  AR,  NF,  FF.  DF.  SF) 

•  f***t**t*»*M*****t**«t*t*»*»M«»9*t»*4**t******»*»****t»***»t»t******»t*»»»tt*»****»»****l» 

FORTRAN  SUBROUTINE  FOR  ORTHONORMAL  POLYNOMIAL  EVALUATION 
*********************************♦**************************♦*♦******♦*******+**************♦ 

The  mode  of  operation  MO  is  0  for  functions,  1  for  first  derivatives,  and  2  for  second 
derivatives.  The  number  of  discrete  arguments  n  is  given  in  argument  NA.  The  variable 
argument  u  is  given  in  argument  AU.  The  recurrence  coefficients  are  given  in  the 
2m  -  2  array  AR.  The  number  of  functions  m  is  given  in  argument  NF.  The  functions 
are  stored  in  the  m- array  FF.  the  first  derivatives  are  stored  in  the  m- array  DF.  and 
the  second  derivatives  are  stored  in  the  m- array  SF. 


ISOMETRIC  REPRESENTATION 


Analysis 


Let  x,  y  be  functions  of  u  such  that  u  is  proportional  to  the  length  of  the  curve 
which  connects  x.  y.  Approximate  values  of  u  are  derived  from  the  perimeter  of  a 
polygon  which  is  inscribed  within  the  curve.  Let  ut  be  the  ith  discrete  value  of  u,  and 
let  xt,  yt  be  the  ith  discrete  values  of  x,  y.  Let  the  values  of  u  be  limited  to  the  range 

-lSuS  +  l  (70) 


The  values  of  u  are  derived  from  the  equation 


2  vW/-  xt)z  +  (yt,t  -  yt)a 
-  *<)*  +  (y»*i  - 


(71) 


Accuracy  of  the  simulation  requires  closely  spaced  data. 

Let  x,  y  be  approximated  as  power  polynomials  in  u  by  the  equations 


n 


*  =  E 

t=  1 


n 


y  =  E  Vi<Pi(u) 

1=1 


(72) 


where  <pt(u)  is  the  ith  Lagrange  polynomial.  The  coefficients  of  the  polynomial 
approximations  of  x,  y  are  the  sums  of  products  of  the  discrete  values  of  x,  y  and  the 
columns  of  the  matrix  of  coefficients  for  the  functions  <fit(u). 

The  polynomial  approximation  for  x.  y  defines  a  curve  for  which  the  metric  l  is 
defined  by  the  equation 


du  '  \\duj  \du) 


(73) 


The  coefficients  of  a  polynomial  approximation  for  l  are  the  sums  of  products  of  the 
discrete  values  of  I  and  the  columns  of  the  matrix  of  coefficients  for  the  functions 
Improved  isometry  is  achieved  if  u  is  replaced  in  accordance  with  the  iteration 


INI 


1  4-  2  :  ‘  - 


dxV 

du) 


m 


dx 

du 


H2) 


du 


du 


(74) 
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Stability  of  the  iteration  requires  closely  spaced  data. 

Programming 

SUBROUTINE  ISMTRS  (NA,  AX.  AY,  SU,  AU) 

t******************************************************************************************** 

FORTRAN  SUBROUTINE  FOR  ISOMETRIC  POLYGONAL  SIMULATION 

v******************************************************************************************** 

The  number  n  of  data  is  given  in  argument  NA.  The  coordinates  x,  are  given  in  the 
n-array  AX,  and  the  coordinates  yi  are  given  in  the  n- array  AY.  The  perimeter  of  the 
polygon  is  stored  in  address  SU.  The  coordinates  ut  are  stored  in  the  n-array  AU. 

SUBROUTINE  LPLNMA  (NA.  AU,  AX,  AY,  CX.  CY,  AC) 

t********************************************************.*************************.*********** 

FORTRAN  SUBROUTINE  FOR  LAGRANGIAN  POLYNOMIAL  APPROXIMATION 

t******************************************************************************************** 

The  number  n  of  data  is  given  in  the  address  NA.  The  arguments  u,,  -,u„  are  given 
in  the  n- array  AU,  the  arguments  x,.--,xn  are  given  in  the  n- array  AX,  and  the 
arguments  are  given  in  the  n-array  AY.  The  coefficients  of  the  polynomial 

for  x  are  stored  in  the  n- array  CX,  and  the  coefficients  of  the  polynomial  for  y  are 
stored  in  the  n-array  CY.  The  coefficients  of  the  Lagrange  polynomials  are  stored  in 
the  (n  +  1)  x  n  array  AC.  A  call  is  made  to  Subroutine  LGRNGX. 

SUBROUTINE  OPlNMA  (NA.  AU,  AX,  AY,  NC,  CX,  CY,  AA,  AR,  AC) 

V******************************************************************************************** 

FORTRAN  SUBROUTINE  FOR  ORTHONORMAL  POLYNOMIAL  APPROXIMATION 

t^ft***!**************,*;*********************************************************************** 

The  number  n  of  data  is  given  in  the  address  NA.  The  arguments  u.,,-  .u^  are  given 
in  the  n-array  AU,  the  arguments  x,,  "-.xn  are  given  in  the  n-array  AX,  and  the 
arguments  y,,  •••,y„  are  given  in  the  n-array  AY.  The  number  m  of  coefficients  is  given 
in  the  address  NC.  The  coefficients  of  the  polynomial  for  x  are  stored  in  the  m-array 
CX,  and  the  coefficients  of  the  polynomial  for  y  are  stored  in  the  m-array  CY.  Values 
of  orthonormal  polynomials  are  stored  in  the  m  x  n  array  AA,  the  recurrence  coefficients 
are  stored  in  the  2m  2  array  AR,  and  the  coefficients  of  the  polynomials  are  stored 
in  the  m  *m  array  AC.  Calls  are  made  to  Subroutine  ORTHOS  and  to  Subroutine  ORTHOX. 

SUBROUTINE  ISMTRL  (NA,  AU,  AL.  CX,  CY.  CL,  AC) 

FORTRAN  SUBROUTINE  FOR  ISOMETRIC  LAGRANGIAN  ITERATION 

The  number  of  arguments  n  is  given  in  the  argument  NA.  The  coordinates  it,,  ■  .  u„ 
are  given  in  the  n  array  Au.  The  values  of  the  l  metric  are  stored  in  the  n  array 
Al.  The  coefficients  for  the  x  -  coordinate  are  given  in  the  n-  array  CX,  and  the  coefficients 
for  the  y  coordinate  are  given  in  the  n-  array  CY.  The  coefficients  of  the  l- metric 
are  stored  in  the  n  array  Cl.  The  matrix  of  Lagrange  polynomials  is  given  in  the 
(n  f  1)  '  n  array  AC.  Improved  coordinates  u,,  -  ,  un  arc  returned  to  the  n  array  AU. 
Calls  are  made  to  Subroutine  MPLNMV. 


SUBROUTINE  ISMTRO  (NA,  AU.  AL,  NC,  CX.  CY.  CL,  AA,  AR,  AC) 

4«**.********************#»*******»********  ******  t***************************  ♦*♦***»*  ***  4***4* 

FORTRAN  SUBROUTINE  FOR  ISOMETRIC  ORTHONORMAL  ITERATION 

t*****************************************************************************************.*** 

The  number  of  arguments  n  is  given  in  the  argument  NA.  The  coordinates  u^.  -,un 
are  given  in  the  n-array  AU.  The  values  of  the  I-metric  are  stored  in  the  n-array 
AL.  The  number  of  coefficients  m  is  given  in  the  argument  NC.  The  coefficients  for  the 
z-coordinate  are  given  in  the  m-array  CX,  and  the  coefficients  for  the  y-coordinate 
are  given  in  the  m-array  CY.  The  coefficients  of  the  I-metric  are  stored  in  the  m-array 
CL.  Values  of  orthonormal  polynomials  are  given  in  the  m  xn  array  AA,  the  recurrence 
coefficients  are  given  in  the  2m  -  2  array  AR,  and  the  matrix  of  coefficients  is  given 
in  the  mxm  array  AC.  Improved  coordinates  ut,  un  are  returned  to  the  n-array 
AU.  Calls  are  made  to  subroutine  MPLNMV. 


LAGRANGE  INTEGRATION 


Analysis 

Let  ult  ••\-un  be  arbitrary  coordinates  and  let/tu,),  ••■,/(un)  be  values  of  an  arbitrary 
function.  Lagrange  interpolation  approximates  the  function  in  accordance  with  the 
equation 


f(u)  =  E  fM^u) 


(75) 


where  <Pt(u)  is  the  ith  Lagrange  polynomial.  Expansion  of  each  Lagrange  polynomial 
in  powers  of  the  argument  and  integration  lead  to  the  integration  multipliers 


m4 


=  J*‘  *(«) 


du 


Then  the  integral  of  an  arbitrary  function  is  given  by  the  equation 

1  n 

/(u)  du  =  E  mt/(wt) 

<=  1 

Expansion  and  integration  is  achieved  by  reference  to  auxiliary  subroutines. 


r.' 


(76) 


(77) 


Programming 

SUBROUTINE  LGRMLT  (NA.  AU.  AM,  CP) 

FORTRAN  SUBROUTINE  FOR  LAGRANGE  INTEGRATION  MULTIPLIERS 

t*************^****^*********************************************************************** 

The  number  of  coordinates  n  is  given  in  argument  NA.  The  coordinates  14,  •■■,un  are 
given  in  the  n-array  AU.  The  integration  multipliers  m,,  mn  are  stored  in  the  n-array 
AM.  The  coefficients  of  the  Lagrange  polynomials  are  stored  in  the  (n  +  1)  ■  n  array 
CP.  Calls  are  made  to  Subroutine  LGRNGX  and  Subroutine  MPLNMV. 
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ORTHONORMAL  INTEGRATION 


Analysis 

Let  u,.  •••,  un  be  arbitrary  coordinates  and  let/(uj),  -,/(uj)  be  values  of  an  arbitrary 
function.  Orthonormal  interpolation  approximates  the  function  in  accordance  with  the 
equation 

f(u)  =  E  E  (7s) 

t=0  i=l 

where  <pt(u)  is  the  fcth  orthonormal  polynomial.  Expansion  of  each  orthornormal 
polynomial  in  powers  of  the  argument  and  integration  lead  to  the  integration  multipliers 

n-l  p+1 

=  E  <Pt(u t)  VtM  du  (79) 

it' 0  J-l 

Then  the  integral  of  an  arbitrary  function  is  given  by  the  equation 

J+l  n 

f(u)  du  =  £  rriiffUi)  (80) 

-l  <=i 

Expansion  and  integration  is  achieved  by  reference  to  auxiliary  subroutines. 


Programming 


SUBROUTINE  OPLMLT  (NA,  AU,  AM,  AA.  AR,  AC) 

***#****i|l**************  +  *  +  *****  +  ****  +  *  +  +  **  +  ****  +  ***-*  +  *+  *  +  +  *+  ****  +  **  +  +  *'**.*Ak**A»*.4.  A  ♦  t.  •***»*<♦* 

FORTRAN  SUBROUTINE  FOR  ORTHONORMAL  POLYNOMIAL  INTEGRATION  MULTIPLIERS 


The  number  of  coordinates  n  is  given  in  argument  NA.  The  coordinates  Uj.-  .u*  are 
given  in  the  n-  array  AU.  The  integration  multipliers  mn  are  stored  in  the  n  -array 

AV.  The  values  of  orthonormal  polynomials  are  stored  in  the  n  x  n  array  AA,  the 
recurrence  coefficients  are  stored  in  the  2n  •  2  array  AR,  and  the  coefficients  of 
expanded  polynomials  are  stored  in  the  n  x  n  array  AC.  Calls  are  made  to  Subroutine 
ORTHOS,  Subroutine  ORTHOX,  and  Subroutine  MPLNMV. 


TRAPEZOIDAL  INTEGRATION 


Analysis 

Let  0,,  •••,  0n  be  equally  spaced  values  of  an  angle  0  such  that  the  ith  value  is  defined 
by  the  equation 


0t  =  (i~D—  (81) 

n 

The  Euler  theorem  and  the  rule  of  summation  for  a  geometric  series  are  used  in  the 
summation  of  the  products  of  the  trigonometric  functions 

cos  fc0  sin  md  (82) 


18 


•A-  ,9*  .  At 


I 

'1  . 


The  summations  over  the  angles  flt  are  given  by  the  equations 

~  „  „  sin(fc  -  m) 271  sin(fc  +  m)Z-n  .  ,  . 

Y.  cos  kdi  cos  m8t  = - - - —  +  — — - —  (i fc  ±  to  *  in)  (83) 

t7i  4  sin(fc  -  to)£  4  sin(fc  +  m)£  '  '  '  '  ' 

n 

cos  fc<?4  sin  to04  =  0  (84) 

i»l 

”  .  .  „  „  sin(fc  -  m)2Tr  sin(fc  +  m)27T  ,  ,  , 

T.  sin  fcfl4  sin  mdt  =  - ; - — - — - —  (fc  ±  m  *  jn)  (85) 

4  sin (k  -  m)l  4  sin(k  +  m)£  v  J  ’  v  ; 


and  by  the  equations 


X)  cosafc04  =  n 


^  cosafc04  =  |n 

i=  l 
n 

J  sinafc04  =  |n 


(fc  =  0)  (86) 


(0  <k<{n)  (87) 


(0  <  fc  ^  |n)  (88) 


The  functions  are  orthogonal  if  fc,  m  satisfy  the  inequality 

0  <  |fc  ±  m|  <  n  (89) 

No  term  with  fc  =  in  can  be  included  in  the  cosine  series  because  cos  fc04  is  zero  for 
all  i  when  fc  =  |n.  An  arbitrary  function  /(0)  is  approximated  by  the  trigonometric 
polynomial 

m<in  m*  |n 

/(B)  =  +  £  oBcosmfl+  £  bmSinme  (90) 

m=l  m=l 

The  coefficients  in  the  approximation  are  recovered  with  the  equations 


“m  =  “  S  /(0 t)  cos  m04 

^  *=t 


=  -  E  /(0t)  sin  to04 

^  4=1 


(m  <  in)  (91) 


(m  5  -J-n)  (92) 


Only  the  constant  Og  survives  integration  through  the  angle  2 n  as  expressed  by  the 
equation 

C  md8^-~if(8 4)  (93) 

Jo  71  i=l 

which  is  the  trapezoidal  rule  of  integration  over  a  complete  cycle. 

Programming 

SUBROUTINE  TPZMIT  (NA,  AO,  AM) 

********************************4*******4**4i**4*******4i*****A*44*4.********.*****«i*****4.*»***t* 

FORTRAN  SUBROUTINE  FOR  TRAPEZOIDAL  INTEGRATION  MULTIPLIERS 

*********************************A******4L***«****4*4A4*4**44***44*»4*4-****»**A*.4»aA4  *aAAA*AA* 

The  number  n  of  coordinates  is  given  in  argument  NA.  The  equally  spaced  coordinates 


The  summations  over  the  angles  04  are  given  by  the  equations 

"  sin(fc  -  m)2n  sin(fc  +  m)2n 

Y  COS  fc04  COS  7710 1  = : — t- - rr  +  

4„!  4  sin(fc  -  tti)-  4  sm(fc  +  m)J 

n 

Y  cos  fc04  sin  77i0t  =  0 

*=i 

sin(fc  -  77i)2tt  sin(fc  +  m)27r 


Y  sin  fc04  sin  m04  = 


4  sin(fc  -  771.)“  4  sin(fc  +  m)% 


and  by  the  equations 


Y  cos 2fc04  =  n 

t=i 
n 

2  cos zkQi  -  |n 

i=  1 

n 

2  sin2fc04  =  |n 

«=i 

The  functions  are  orthogonal  if  k,  m  satisfy  the  inequality 

0  <  ]fc  ±  77l|  <  71 

No  term  with  fc  -  jn  can  be  included  in  the  cosine  series  because  cos  fc04  is  zero  for 
all  i  when  k  =  ^ n .  An  arbitrary  function  /(0)  is  approximated  by  the  trigonometric 
polynomial 


(fc  ±  77i  ^  jn)  (83) 
(84) 

(fc  ±  m  *  jn)  (85) 

(fc  =  0)  (86) 
(0  <  fc  <  |t i)  (87) 
(0  <  fc  g  |ti)  (88) 

(89) 


mCgU 

/(0)  =  |a0  +  £  am  cos  7710+  X)  6m  sin  77i0 

m=l  m=*  1 


(90) 


The  coefficients  in  the  approximation  arc  recovered  with  the  equations 

2  n 

=  -  Y  COS  77l0t  (tti  <  |n)  (91) 


bm  =  -  £  f(0i)  sin  7ti04 


71 


( 771  S  |n)  (92) 


Only  the  constant  a0  survives  integration  through  the  angle  277  as  expressed  by  the 
equation 


f 


m  dg  -=  27T  Y  Adi) 


(93) 


which  is  the  trapezoidal  rule  of  integration  over  a  complete  cycle. 
Programming 

Subroutine  tpzmlt  (na,  ao,  am) 

**************************  +  **************  +  *♦****%****.*****.  +  ***♦********♦*  +  *,*.***** 
FORTRAN  SUBROUTINE  FOR  TRAPEZOIDAL  INTEGRATION  MULTIPLIERS 


******.*** 


The  number  ti  of  coordinates  is  given  in  argument  NA.  The  equally  spaced  coordinates 
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0t.  0n  are  stored  in  the  n-array  AQ.  The  integration  multipliers  ml,-  -,mn  for 
trapezoidal  integration  are  stored  in  the  n-array  AM. 


CHEBYSHEV  INTEGRATION 

Analysis 

Let  Bl,-,8n  be  equally  spaced  values  of  the  angle  0  such  that  the  ith  value  is 
defined  by  the  equation 


n 


(94) 


The  Euler  theorem  and  the  rule  of  summation  for  a  geometric  series  are  used  in  the 
summation  of  the  products  of  the  trigonometric  functions 


cos  kB  cos  m3 

The  summations  over  the  angles  0t  are  given  by  the  equations 


(95) 


"  „  sin(fc  -  m)7T  sin(fc  +  m)iT  ,  ,  , 

V  cos  kBi  cos  mOi  = - - - —  +  - - — - —  (k  ±m*  2 in)  (96) 

t.j  ‘  1  4  sin(fc  -  m)/n  4  sin(A:  f  m)£  v  J  ’ 


and  by  the  equations 


Yj  cos zk6i  =-  n 

t=i 

(k 

0)  (97) 

n 

V  coszfc0t  =  |n 

•i=l 

The  functions  are  orthogonal  if  k,  m  satisfy  the  inequality 

(0  <  k  < 

n)  (98) 

0  <  \k  ±  m!  <  2n 

(99) 

An  arbitray  function  /(0)  is  approximated  by  the  trigonometric  polynomial 

n- 1 

/(0)  =  |a0  +  am  cos  77x0 

m~  i 

(100) 

The  coefficients  in  the  approximation  are  recovered  with  the  equations 

2  n 

a-m  -  ~  I!  /(6 i)  cos  m0t 


(101) 


Only  the  constant  a0  survives  integration  through  the  angle  n  as  expressed  by  the 
equation 


f 


m  do  -V  /(0t) 


which  is  the  trapezoidal  rule  of  integration  over  a  half  cycle. 

Let  the  argument  x  be  defined  in  terms  of  the  angle  0  by  the  equation 

x  -  -  cos  0 


(102) 


(103) 
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**  r.n-y&fgf$9QX'-4t) 


Integration  with  respect  to  x  is  expressed  by  the  equation 

T 

/(0)  sin  0  d8 


J+ 1  f*ir 

f(x)  dx  =  J 


(104) 


Term  by  term  integration  of  the  trigonometric  polynomial  is  achieved  with  the  aid  of 
the  indefinite  integral 


I 


cos  k6  sin  0  d6 


cos (fc  -  1)0  cos(A:  +  1)0 

2(k  l)  Z(kVlT 

Evaluation  at  the  limits  of  integration  leads  to  the  equations 


f 

Jo 

f 

Jo 


cos  kd  sin  d  d6  -  0 


cos  kd  sin  8  d8  =  - 


jo  (*2-l) 

Integration  multipliers  are  defined  by  the  equation 

2  4  cos  2k6i 

mi  ~  n  n  Jr,  Akz  -  1 

Then  integration  with  respect  to  x  is  evaluated  by  the  equation 


f 


f(x)  dx  =  £  rr^fixj 

1  i~  1 


(105) 

(odd  k)  (106) 
(even  k)  (107) 

(108) 

(109) 


where  the  mt  are  the  integration  multipliers  for  Chebyshev  integration. 
Programming 

SUBROUTINE  CHVMLT  (NA,  AU,  AM) 

FORTRAN  SUBROUTINE  FOR  CHEBYSHEV  INTEGRATION  MULTIPLIERS 

*************4****************** ***444**44 4 *44  4  4  »  *44*44  4.  *  »  A  »  *.  4  4  *  4  4  M  *  *  *14  M«  i  t ‘4  *  *  *  *  *  *  4*44*4* 

The  number  n  of  coordinates  is  given  in  argument  NA.  The  cosines  u,,  •  ■  ,un  of  equally 
spaced  angles  are  stored  in  the  n-  array  AU.  The  integration  multipliers  mt,  •,mn  for 
Chebyshev  integration  are  stored  in  the  n-array  AM. 


HIGH-ACCURACY  QUADRATURE 

Analysis 

Christoffel-Darboux  Identity 

Let  an  arbitrary  function  f(x)  be  given  at  the  abscissae  r,,  -,xn.  Let  the 
polynomial  p(x)  be  equal  to  the  values  of  f(x)  at  the  abscissae  as  expressed  by  the 
equation 

p(*t)  =/(*«)  ("0) 
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Let  the  function  F(x)  be  defined  by  the  equation 

F(x)  =  (x  z,)-(z-zn)  (111) 

The  function  F(x)  is  a  power  polynomial  of  the  nth  degree  and  is  zero  at  each  abscissa. 
Let  the  function  q(x)  be  a  power  series  which  satisfies  the  equation 

/(z)  =  p(x)  +  q(x)F(x)  (112) 


If  the  function /(z)  itself  were  expressed  by  a  power  series,  then  the  coefficients  of  q{x) 
could  be  determined  one  by  one  by  an  algorithm  which  starts  at  the  highest  power 
of  z  and  compares  terms  with  the  same  power  of  z.  The  abscissae  are  so  selected  as 
to  satisfy  the  equation 


f 


q(x)F{x)w(x)  dx  =  0 


(113) 


where  ui(z)  is  a  weighting  function,  and  q(x)  is  any  polynomial  of  the  (n  -  l)th  degree. 
Insofar  as  f(x)  can  be  approximated  by  a  polynomial  of  the  (2n  -  l)th  degree,  the 
integral  of  /(z)  is  given  accurately  by  the  equation 


^6 

f(x)w(x)  dx  =  p(x) w(x)  dx 

J  a  v  a 


(114) 


The  function  q(z)  can  be  expressed  as  a  linear  combination  of  orthonormal  polynomials 
of  which  the  nth  polynomial  is  the  function  F(x). 

Let  (p0(z),  $Bn_,(z)  be  n  polynomials  of  progressively  increasing  degree.  The 

polynomials  are  orthonormal  if  they  satisfy  the  condition 

f  <PkM<fim(x)w(x)  dx  =  ^  km  (115) 

J  a 


where  6km  is  zero  if  fc  *  m  and  is  unity  if  fc  =  m.  Any  polynomial  of  degree  m  can  be 
expanded  in  terms  of  the  orthonormal  polynomials  of  degrees  0  through  m.  The 
coefficients  in  the  expansion  are  given  by  the  solution  of  a  sequence  of  equations 

which  start  at  the  highest  degree  and  compare  terms  of  the  same  degree.  Let  <pm(x) 

be  expressed  by  the  equation 

m-1 

<Pm(x)  =  7mx<pm-i(x)  E  7v<Pv(x)  (116) 

u=0 

where  the  yu  are  constant  coefficients.  Multiplication  throughout  by  <pk(x)  and  integration 
isolates  yk  in  accordance  with  the  equation 

pb  pb 

J  ‘Pt(x)vm(x)w(x)  dx  =  ym  J  x<pk(x)?m_i(x)ui(x)  dx  +  7*  =  0  (117) 

The  polynomial  x<pk(x)  may  be  expressed  by  the  equation 

fc+  i 

x<pk(x)  =  E  e„p„(x)  (118) 

v=0 


where  the  e„  are  constant  coefficients.  All  terms  of  the  summation  are  orthogonal  to 
<Pm  - i(x)  an<l  the  coefficient  yk  is  zero  for  all  fc  S  m  -  3.  The  orthonormal  polynomials 
can  be  generated  by  a  three  term  recurrence  equation. 
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If  is  zero  and  <p0(x)  is  constant,  then  the  three-term  recurrence  equation  is 

<PmW  =  — —  ~  Pm-lVm-dx)  ~  am-2<fim-Z(x)\  (119) 

where  the  constants  am.i  and  0m-,  are  defined  by  the  equations 


J  a 

—r 


XVm-l  (x)<pm(x)w(x)  dx 


X(fim-Ax)<pm  ^x)ui(x)  dx 


Multiplication  of  the  three-term  recurrence  equation  throughout  by  cpm._  t(y), 
interchange  of  x  and  y,  and  subtraction  gives  the  equation 

(x  -  y)<Pm-l(x)Vm-l(y)  ~  «„■  !  lvOm(*Vm-l(?/)  ~  <Pm- 1  (x)<fim(y)  j 

-  am-Z  \<Pm-l(x)<Pm  z(y)  ~  <fm. 1  (y )  !  (122) 

Summation  with  respect  to  m  gives  the  Christoffel-Darboux  identity 

n  - 1 

«n-l  1  <Pn(x)<fin-l(y)  -  ¥>„-i(*)^n(l/)l  =  (*  ~  y)  T.  <Pt(x)<pt(y)  (123) 


Rearrangement  and  integration  gives  the  equation 

f  *,(*)  \  ±n(x)9n-i(y)-9n->(x)Vn(y)  )  ^  =  fjM  (124) 

Ja  (  x  y  )  a„  , 

which  is  fundamental  to  Lagrange  integration  when  y  is  set  equal  to  an  abscissa  a:*. 

Chcbyshev  Quadrature 

The  Chebyshev  polynomial  Tn(x)  is  defined  by  the  equation 

Tn(x)  -  cos(n  cos  ‘x)  (125) 

The  Chcbyshev  polynomial  for  0  v  n  is  given  by  the  equation 


T  lx)  2"  •*  V  (  Dmn(n  -  m)!»"  ^ 
2Zm(n  -  m)(n  -  2 m)!m! 


The  functions  of  lowest  order  are  given  by  the  equations 


T0(x)  1  T,(x)  =  x 

and  the  functions  satisfy  the  recurrence  equation 

TJx)  «  2xTn  .t(x)  -  Tn  z(x) 

The  first  derivatives  of  lowest  order  arc  given  by  the  equations 

To(x)  0  T[(x)  1 

and  the  first  derivatives  satisfy  the  recurrence  equation 


(0  <  n)  (126) 


T'n(x)  -  nTn  ,(x)  -r  --  -*7V,(x) 


The  second  derivatives  of  lowest  order  are  given  by  the  equations 
To(x)  =  0  T['(x)  -  0 

and  the  second  derivatives  satisfy  the  recurrence  equation 

K(x)  n*-  ■  7'„  ,(x)  *  nTxr;  ,(x) 
n  1  n  -  1 

The  Chebyshev  functions  satisfy  the  orthogonality  equations 


(131) 


(132) 


r 

j- 

Tk(x)Tm(x) 

'  V'l  x2  dl  ° 

(A:  *  to)  (133) 

i: 

W  dx  J’ 

Vl  -  x2  \  >rr 

(n  0) 

(134) 

(O  n) 

The  limits  of  integration  and  the  weighting  function  are  given  by  the  equations 

,  ,  1 

a  ■■■  -  1  te(x)  -  -----  6  ;  1 

v  1  x2 

The  orthonormal  functions  and  the  Chebyshev  functions  are  related  by  the  equation 


(135) 


?n(*)  \J  n  Tn(x) 


(136) 


A  comparison  between  the  three  term  recurrence  equations  leads  to  the  equations 


Qn  1  a 

The  function  F(x)  is  given  by  the  equation 


,  0 


r(x)  21  ,  Tn(x) 


and  its  derivative  is  given  by  the  equation 


r(x)  -  -i  jn(x) 

The  abscissae  arc  given  by  the  equation 

,  .  i-TT 

x4  cos(r  i) 

n 

The  Lagrangian  integration  multipliers  arc  given  by  the  equation 

r)  dx 


m, 


F(x) 

J  ,  F'M 


i)  (x  Xi) 

Substitutions  in  the  ChristofTcl  Darboux  identity  lead  to  the  equation 

rr 

m‘  Tn  ,(xi)r;(xi) 


(137) 


(138) 


(139) 


(140) 


(141) 


(142) 


2-: 


A  comparison  between  the  three  term  recurrence  equations  leads  to  the  equations 


v(2n  -  l)(2n  -f  1) 


fin-l  =  0 


The  function  F(x)  is  given  by  the  equation 


and  its  derivative  is  given  by  the  equation 

The  Lagrangian  integration  multipliers  are  given  by  the  equation 


r  f(x) 

**  J-,  *"(*,) 


dx 


«)  (*  -  *<) 

Substitutions  in  the  Christoffel  -Darboux  identity  lead  to  the  equation 

2 


m,  - 


nPn-  A*i)P'n(Xi) 


The  recurrence  equation 


(1  -  xz)P’n(x)  -  n(Pn..,(x)  xPn(x)) 


leads  further  to  the  equation 


m< "  (T:  x/)|p;(i<)l* 

The  integral  of  an  arbitrary  function  is  given  by  the  equation 


f  f(l)dX  -  Y. 

J-i  t=l 


where  the  mt  arc  the  multipliers  for  Gaussian  integration. 

Let  0(  be  an  angle  which  is  defined  by  the  equation 

i  -  ! 

0 i  -  " 

n  l-  g 

Then  an  initial  approximation  for  x(  is  given  by  the  equation 

x<  ~  ±  cos  9i 

The  approximation  for  xi  is  refined  by  Newton  Raphson  iteration  until 

PnM  =  0 

to  within  rounding  error.  The  Legendre  function  and  its  derivatives  are 
recurrence  equations. 


(156) 

(157) 

(158) 

(159) 

(160) 

(161) 

(162) 

(163) 

(164) 

(165) 

(166) 
computed  by 
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Programming 

SUBROUTINE  GSSMLT  (NA,  AU,  AW) 

*******iii*********«**«********  +  ****«*»»*********  +  **************4.****»****t********t*********** 

FORTRAN  SUBROUTINE  FOR  GAUSS  INTEGRATION  MULTIPLIERS 

******************************  ********  A****************************************************** 

The  number  n  of  coordinates  is  given  in  argument  NA.  The  roots  u,.  ■■■,un  of  the  nth 
degree  Legendre  polynomial  are  stored  in  the  n  array  AU.  The  integration  multipliers 
mlt  -  ,mn  are  stored  in  the  n  -array  AM. 


COMPLEX  POWER  POLYNOMIAL  EVALUATION 

Analysis 

Let  f(z)  be  a  polynomial  of  (n  -  l)th  degree  in  the  argument  z.  Let  the  polynomial 
be  expressed  in  terms  of  its  argument  by  the  equation 


f(z)  =  £  atzk 


k-0 


(167) 


where  the  coefficients  ak  are  complex.  The  coefficients  are  arranged  in  ascending  order 
in  an  array  with  real  parts  in  alternate  addresses  and  with  imaginary  parts  in  the 
next  higher  addresses.  The  first  derivative  of  the  polynomial  is  given  by  the  equation 


df(z) 

dz 


n-1 


£  katzk  1 

k-0 


(168) 


and  the  second  derivative  of  the  polynomial  is  given  by  the  equation 

d~-T-  =  £l*(*- Da**4'*  (169) 

dz  *=o 

Whether  the  operation  is  evaluation,  first  differentiation,  or  second  differentiation  is 
determined  by  a  mode  of -operation  parameter.  The  polynomials  are  evaluated  by  the 
nested  method  of  polynomial  evaluation. 


Programming 

SuBROd'Nr  CPWRV  (WO,  A,:,  NC.  AC,  FV) 

***********•*****♦*•********♦***♦*♦**»♦#♦**♦*♦+**♦**♦♦*♦******♦#***♦***♦*******♦♦*********%•.* 
FORTRAN  SUBROUTINE  FOR  COMPLEX  POWER  POLYNOMIAL  EVALUATION 

The  mode  of  operation  is  evaluation  if  MO  --  0,  first  differentiation  if  MO  =  1,  and  second 
differentiation  if  MO  =  2.  The  argument  z  is  given  in  the  2  array  AZ.  The  number  of 
coefficients  n  is  given  in  the  address  NC.  The  coefficients  a0,  ■••,aT,_l  are  given  in  the 
2 n  array  AC.  The  value  of  the  function  is  stored  in  the  2  -array  FV. 


k 
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COMPLEX  FOURIER  SERIES  EXPANSION 


Analysis 

Let  0l.  -.0n  be  a  set  of  angles  such  that 

0*  =  (*-!)—  (HO) 

n 

Each  angle  is  the  midpoint  of  one  of  n  intervals  into  which  2n  is  divided.  Let  a 
function  f(6)  be  expressed  in  terms  of  its  argument  0  by  the  direct  transform 

(171) 

k  =  0 

where  the  coefficients  A*  are  complex.  The  coefficients  are  recovered  from  values  of 
the  function  by  the  inverse  transform 

Am  -  '  V  /(0t)c"‘m,,‘  (172) 

n  t 

The  weighting  factors  for  the  transform  are  stored  in  an  n  '  2m  matrix  of  coefficients. 
The  weighting  factors  are  the  same  for  all  values  of  m(k  -  lz)  which  differ  by  n.  It  is 
sufficient  to  compute  the  2n  pairs  of  trigonometric  functions  for  which  k  1  and  k  -  | 
are  less  than  n.  The  trigonometric  functions  are  computed  and  are  stored  in  the  first 
two  columns  of  the  matrix.  The  remainder  of  the  matrix  is  filled  by  transfers  from 
the  first  two  columns.  Then  the  first  column  is  replaced  by  constants. 

The  coefficients  for  a  particular  set  of  values  of  the  function  are  obtained  from  a 
vector -matrix  multiplication.  From  the  Fourier  transform  for  one  set  of  angles  may 
be  derived  the  Fourier  transform  for  another  set  of  angles  If  the  angles  are  displaced 
by  a  constant  c,  then  the  coefficients  of  the  Fourier  transform  arc  multiplied  by  powers 
of  e“ 

Programming 

SUBROUTINE  CFOUPX  (NA.  NC,  AC) 

************ ^*******#*****************#*****A***.»*****»***#**»*******#****»  **•****►*»**♦•♦•** 

FORTRAN  SUBROUTINE  FOR  COMPLEX  FOURIER  SERIES  EXPANSION 
*******♦*************%»*******♦*****•*•****♦**************#*******♦******♦•**♦****•**•****•»• 

The  number  n  of  data  is  given  in  argument  NA,  and  the  number  m  of  weighting  factors 
is  given  in  argument  NC.  The  matrix  of  weighting  factors  is  stored  in  the  n  ■  2m 
array  AC. 


COMPLEX  FOURIER  SERIES  EVALUATION 


Analysis 

Let  /(0)  be  a  polynomial  of  (n  -  l)th  order  with  angular  argument  0.  Let  the 
polynomial  be  expressed  in  terms  of  its  argument  by  the  equation 

m  =  V  ^eime  (173) 

TO"  0 

where  the  coefficients  Am  are  complex.  The  coefficients  are  arranged  in  ascending 


28 


T 


order  in  an  array  with  real  parts  in  alternate  addresses  and  with  imaginary  parts  in 
the  next  higher  addresses.  The  first  derivative  of  the  polynomial  is  given  by  the 
equation 


dm 

de 


n-1 

S.  im6 

xrnAme 

m~0 


(174) 


and  the  second  derivative  of  the  polynomial  is  given  by  the  equation 


dzf{6) 

dBz 


-  T.  ™24me‘' 


(175) 


Whether  the  operation  is  evaluation,  first  differentiation,  or  second  differentiation  is 
determined  by  a  mode  of-operation  parameter.  The  polynomials  are  evaluated  by  the 
nested  method  of  polynomial  evaluation. 


Programming 

SJ3R01TINE  CFOURV  (MO.  AQ,  NA,  AA,  FA/) 

FORTRAN  SUBROUTINE  FOR  COMPLEX  FOURIER  SERIES  EVALUATION 


The  mode  of  operation  is  evaluation  if  MO  =  0,  first  differentiation  if  MO  -  1,  and  second 
differentiation  if  MO  =  2.  The  argument  0  is  given  in  the  address  AQ.  The  number  of 
coefficients  n  is  given  in  the  address  NA.  The  coefficients  A0.  ■■■ ,An  i  are  given  in  the 
2 n  array  AA.  The  value  of  the  function  is  stored  in  the  2- array  FV. 


FAST  FOURIER  TRANSFORM 


Analysis 

In  the  conventional  discrete  transform  a  cycle  is  divided  into  N  intervals.  A  set  of 
discrete  arguments  60,  is  defined  by  the  equation 

2tt 

ek-k—  (176) 

A  function  /(0)  is  expressed  in  terms  of  the  argument  6  by  the  series  in  the  equation 

m  =  j:  Ameimt  (177) 

m=0 

The  values  of  the  function  at  the  discrete  arguments  are  given  by  the  discrete  Fourier 
transform  in  the  equation 

/(«*)  =  *£  (178) 

m-0 

The  formula  for  the  summation  of  geometric  series  leads  to  an  orthogonality  relation 
as  expressed  by  the  equation 

^  -  N6nm  (179) 

*=0 
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where  Snm  is  zero  if  m  *  n  and  is  unity  if  m  =  n.  The  coefficients  of  the  Fourier  series 
are  isolated  by  applications  of  the  orthogonality  relation.  The  coefficients  are  given 
by  the  equation 

1  N1  -km—i 

Am  =  ~Zetm»if(8k)  (1B0) 

N  t=o 

In  the  summation  with  respect  to  k,  there  is  redundancy  when  N  is  factorable. 

If  N  =  2n,  then  the  indices  k,  m  are  expressed  as  binary  numbers  in  accordance 
with  the  equations 

fc  =  fc„.12nI+-  +  fc0l  (181) 

m  =  mn_12n_1  +  +  to0  ■  1  (182) 


where  each  binary  bit  has  the  values  0,  1.  Single  summation  with  respect  to  k  is 
replaced  by  multiple  summation  with  respect  to  fcn_,  .■■■,k0  and  functions  of  k  or  m 
are  replaced  by  functions  of  .  k0  or  mn.l,  Thus  the  Fourier  transform  is 

expressed  by  the  equation 


-.mo)  =  ^  £  •  £  /(fc„-,,  -.fc0)^<‘n“l2 

■N  t  -0  k  =0 

n- 1  0 

where  the  weight  w  is  defined  by  the  equation 

xu  —  e  K 


+*0>(mn-l2n"  1+'+m0) 


(183) 


(184) 


The  weight  w  satisfies  the  identities 

•ui*  =  1  to2*  =  -  1  (185) 

There  is  therefore  a  repetition  of  factors  whenever  km  is  incremented  by  a  multiple 
of  ±<V. 

In  the  Cooley-Tukey  version  of  the  FFT,  the  binary  expression  for  k  is  resolved  into 
its  individual  terms.  In  the  product 

(mn^2n-'  +  -  +  m0)kv-l2v~1  (186) 

all  terms  reduce  w  to  unity  except  those  in  the  product 

(ran_l/2n~l/  +  •••  +  m0)kv-l2v>  (1B7) 


The  weighting  factor  in  the  Fourier  transform  is  given  by  the  equation 
wkm  =  xu(m°)kn  i2n~‘+  -+<m«-i2n  l+  •+mo)to 


(188) 


The  summations  with  respect  to  ifen_ t ,  •••,  k0  are  nested.  The  summation  with  respect 
to  kn_v  is  associated  with  m,,., 2v~l  +  •••  +  m0.  If  the  transformations  of  data  are  made 
in  place,  then  new  arrays  of  data  are  created  with  the  aid  of  the  recurrence  equation 


''*0)  ~  £ 


•.t0)ie 


(189) 


where  the  order  of  the  bits  Tnn_1,  -,m0  becomes  reversed.  The  original  set  of  data  is 
transformed  into  new  sets  of  data  until  the  final  set  is  the  set  of  coefficients.  Initially, 
summation  is  for  kn.i  =  0,  1  with  all  values  of  fcn-2,  •  ,  k0  and  with  m ^  -  0.  1.  Finally, 
summation  is  for  k0  -  0,  1  with  all  values  of  m0,  .  mn_t  In  each  cycle  of  transformation 
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the  number  of  sets  of  data  is  doubled  while  the  number  of  data  per  set  is  halved. 

In  the  Sande-Tukey  version  of  the  FFT,  the  binary  expression  for  m  is  resolved  into 
its  individual  terms.  In  the  product 

(fcn_12n_1  +  +  fc0)m„_121'~1  (190) 

all  terms  reduce  w  to  unity  except  those  in  the  product 

(fcn_1/2n_1'  +  +  fc0)ml/_l2,'~1  (191) 

The  weighting  factor  in  the  Fourier  transform  is  given  by  the  equation 

.^*m  =  lu<*o>mn-i2n“1+ -+<*«-i*n'1+-+*o>mo  (j92) 

The  summations  with  respect  to  fcn_(.  k0  are  nested  with  the  summation  with  respect 
to  kn  v  associated  with  m„_,.  If  the  transformations  are  made  in  place,  then  new 
arrays  of  data  are  created  with  the  aid  of  the  recurrence  equation 

while  the  order  of  the  bits  ■  ,m0  becomes  reversed.  The  original  set  of  data  is 

transformed  into  new  sets  of  data  until  the  final  set  is  the  set  of  coefficients.  Initially, 
summation  is  for  kn  ,  -  0,  1  with  all  values  of  fcn_2,  ••• ,  k0  and  with  m0  -  0,  1.  Finally, 
summation  is  for  k0  -  0,  1  with  all  values  of  m0,  In  each  cycle  of  tranformation 

the  number  of  sets  of  data  is  doubled  and  the  number  of  data  per  set  is  halved. 

The  difference  between  the  two  algorithms  is  in  the  powers  of  w  which  are  associated 
with  summations  with  respect  to  the  binary  bits.  Both  algorithms  require  nlog2n  cycles 
of  transformation,  and  both  are  unscrambled  in  accordance  with  bit  reversal  if  they 
are  excecuted  in  place. 

Programming 

SJ-A'^Oy  *.*  C-:  :  •  (VO,  Ji,  AA) 

FORTRAN  St'BROl  TINE  FOR  COMPLEX  FAST  FOURIER  TRANSFORM 
••*•••••••••■•••«•««««•««•**+»****»••#*»#***+******************************************♦***** 

The  mode  of  operation  is  given  in  argument  MO.  The  transform  is  from  data  to 
coefficients  if  V  '  1.  and  from  coefficients  to  data  if  MO  =  +1.  The  log2A'  is  given  in 

argument  '.  The  data  or  coefficients  are  given  or  stored  in  the  2N  array  AA. 


FORTRAN  SUBROUTINE  FOR  BIT  REVERSAL 


multiples  of  (1,  N)2n  for  a  quadrant  are  stored  in  the  | N  +  2  array  WK. 


SUBROUTINE  MXFFT  (MO,  LN.  IR,  WK,  IA,  AA) 

«**«^******4^**%***t^*******************i^************»**l|^*4^*************'t******t********* 

FORTRAN  SUBROUTINE  FOR  FOURIER  TRANSFORM 

****4^********-**********%**********4^*4^************************************************t**** 

The  mode  of  operation  is  given  in  argument  MO.  The  transform  is  from  data  to 
coefficients  if  MO  -  -  1,  and  from  coefficients  to  data  if  MO  +1.  The  log2Af  is  given  in 
argument  LN.  The  bit  reverse  indices  are  given  in  the  \N  array  IR,  and  the  trigonometric 
functions  are  given  in  the  \N  +  2  array  WK.  The  interval  between  the  addresses  of  data 
or  coefficients  is  given  in  argument  IA.  The  data  or  coefficients  are  given  or  stored  in 
the  2.V  array  AA. 


ROOTS  OF  POLYNOMIALS 

The  solution  of  the  quadratic  equation  was  known  to  the  Arabs  in  the  ninth  century. 
The  solution  of  the  eubic  equation  can  be  credited  to  Tartaglia  in  1530  His  solution 
was  published  by  Cardan  in  1545.  The  solution  of  the  quartic  equation  can  be  credited 
to  Ferrari  who  was  a  pupil  of  Cardan's.  He  resolved  the  quartic  polynomial  into  a  pair 
of  quadratic  factors.  Another  solution  of  the  quartic  equation  was  given  by  Euler  in 
1770  He  assumed  that  the  roots  could  be  expressed  by  the  sums  or  differences  of 
three  radicals  A  correlation  of  the  methods  of  solution  with  their  discoverers  is  to 
be  found  in  a  book  by  Burnside  and  Panton14. 

Solutions  of  the  quartic  equation  can  be  classified  as  two  radical15-16  or 
three  radical17.  The  two  methods  must  be  equivalent  analytically,  but  there  is  a 
difference  in  their  accuracy  of  computation. 

In  the  Ferrari  method,  the  roots  of  the  quartic  equation  are  computed  with  the 
most  positive  real  root  of  a  resolvent  cubic  equation,  whereas  in  the  Euler  method, 
the  roots  of  the  quartic  equation  are  computed  with  all  three  roots  of  the  resolvent 
cubic.  In  a  special  case  of  closely  spaced  roots,  rounding  error  might  upset  the  choice 
of  root  for  the  Ferrari  method,  but  merely  would  modify  the  three  roots  for  the  Euler 
method 

Subroutines  have  been  available  in  the  CDC  6700  library  for  the  solution  of  quadratic, 
cubic,  and  quartic  equations.  However,  these  subroutines  do  not  give  roots  in  a  format 
which  is  useful  for  further  complex  arithmetic,  and  they  do  not  give  roots  with  the 
full  accuracy  of  the  computer.  New  and  improved  subroutines  have  been  prepared  in 
the  present  investigation  The  new  subroutines  force  the  product  of  the  roots  to  be 
equal  to  the  constant  in  order  to  reduce  rounding  error  in  the  smallest  root  Previous 
subroutines  for  the  quartic  equation  have  used  the  Ferrari  method,  whereas  there  are 
new  subroutines  now  for  both  the  Ferrari  method  and  the  Euler  method.  In  the  previous 
version  of  the  Ferrari  method,  a  choice  between  alternate  formulae  was  made  only  to 
avoid  division  by  zero,  whereas  in  the  present  version  of  the  Ferrari  method,  the 
choice  is  made  to  achieve  optimum  accuracy. 

The  function  name  SQRi(')  has  been  preempted  in  FORTRAN  for  the  square  root  of 
-.  The  function  name  CBR'(\)  is  used  on  the  UN'IYAC  1108  computer  for  the  cube  root 
of  -  A  new  function  routine  for  the  cube  root  has  been  prepared  for  the  CDC  6700 
computer.  It  obviates  the  inefficiency  of  exponentiation 
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COMPLEX  ROOTS  OF  A  QUADRATIC  EQUATION 


k  . 


Analysis 

Let  a  quadratic  equation  have  real  coefficients  as  expressed  by  the  equation 

azz  +  bz  +  c  =  0 

The  term  in  z  is  eliminated  by  the  substitution 


Then  f  is  a  solution  of  the  equation 


atz  -  —  +  c  =  0 
4a 

and  z  is  given  by  the  well  known  quadratic  formula 

-  6  ±  v62  -  4  a  c 


or  by  the  equivalent  formula 


2c_  __ 

6  -  v'6*  4ac 


which  gives  better  accuracy  when  c  is  small  and  the  sign  of  -6  is  opposite  to  the  sign 
of  the  radical. 

If  the  discriminant  satisfies  the  inequality 

bz  -  4ac  S  0  (199) 

the  roots  are  real,  but  if  the  discriminant  satisfies  the  inequality 

bz  -AachO  (20C) 

the  roots  are  complex  conjugate. 

The  discriminant  is  exactly  zero  only  for  integer  values  of  the  coefficients.  Otherwise 
there  is  a  residual  error  e  either  from  inherent  error  in  the  coefficients  or  from 
rounding  error  in  the  evaluation  of  the  discriminant.  In  either  case  there  remains  an 
error  el  2  in  the  square  root  of  the  discriminant.  It  is  not  possible  to  distinguish  two 
closely  spaced  roots  from  one  double  root  unless  they  have  a  discriminant  larger 
than  e. 

Programming 

SUSROuU’it:  QDCFT  (AC,  RZ) 

*  *  *  *  *  *  *  *  *  *  *  ******  **************************************************************************** 
FORTRAN  SUBROUTINE  FOR  ROOTS  OF  QUADRATIC  POLYNOMIAL 

**»**»********************«********************************************************4********4 

The  real  coefficients  of  the  polynomial  are  given  in  array  AC  in  ascending  order.  The 
roots  of  the  polynomial  are  computed  with  the  quadratic  formula.  The  real  and 
imaginary  parts  of  the  roots  are  stored  in  alternate  addresses  of  array  RZ. 
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COMPLEX  ROOTS  OF  A  CUBIC  EQUATION 

Analysis 


Let  a  cubic  equation  have  real  coefficients  as  expressed  by  the  equation 


z3  +  pzz  +  qz  +  r  =  0 

(201) 

The  term  in  z2  is  eliminated  by  the  substitution 

z  =  t-?- 
3 

(202) 

Then  t  is  a  solution  of  the  equation 

t3  +  at  +  b  =  0 

(203) 

where  the  coefficients  are  defined  by  the  equations 

a  =  |(3 q  -  pz) 

(204) 

b  =  ^(2p3  -  9 pq  +  27r) 

(205) 

In  Tartaglia's  method,  t  is  given  by  the  equation 

t  =  A  +  B 

(206) 

where  the  constants  A,B  are  defined  by  the  equations 

i 

L  2  \  4  27  J 

3 

(207) 

i 

B  =  [  -  —  -  \j—  +  ~  1 

[  2  V  4  27  J 

13 

(208) 

Among  the  three  complex  cube  roots  for  each  parameter  only  those  pairs  are  selected 
for  which 


AB  -  - 


a 

3 


If  the  discriminant  satisfies  the  inequality 


6* 

4 


+ 


-SO 


27 


(209) 


(210) 


there  are  one  real  root  and  two  complex  conjugate  roots  If  the  discriminant  satisfies 
the  inequa  ity 


(211) 


there  are  three  real  roots.  The  constants  A,  B  are  given  by  the  equations 

® 

\  3  ' 

J-  “  -  > 

\  3 


(212) 

(213) 
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where  0  is  defined  by  the  equation 


<t> 


(214) 


For  every  pair  of  constants  A,  B  the  two  other  pairs  are  obtained  through  multiplication 
by  the  factor 


Since  .4  and  B  are  complex  conjugate,  their  sum  is  real. 

Improved  accuracy  is  achieved  when  the  root  of  smallest  absolute  magnitude  is 
replaced  by  the  quotient  of  the  constant  with  reversed  sign  and  the  product  of  the 
two  largest  roots. 

Programming 

SUBROUTINE  CBCRT  (AC.  RZ) 

FORTRAN  SUBROUTINE  FOR  ROOTS  OF  CUBIC  POLYNOMIAL 


The  real  coefficients  of  the  polynomial  are  given  in  array  AC  in  ascending  order.  The 
roots  of  the  polynomial  are  computed  with  the  Tartaglia  method.  References  are  made 
to  function  routine  CBRT.  The  real  and  imaginary  parts  of  the  roots  are  stored  in 
alternate  addresses  of  array  RZ. 


COMPLEX  ROOTS  OF  A  QUARTIC  EQUATION 

Analysis 

Let  a  quartic  equation  have  real  coefficients  as  expressed  by  the  equation 

z*  +  pz3  qzz  +  rz  s  =  0  (216) 

In  Ferrari's  method,  the  quartic  polynomial  is  factored  into  the  product  of  two  quadratic 
polynomials.  The  solutions  of  the  quartic  equation  are  solutions  of  the  quadratic 
equations 

22  +  kpz  +  it  ±  {az  +  b)  =  0  (217) 

Multiplication  of  the  quadratic  factors  and  comparison  of  coefficients  shows  that  the 
coefficients  are  related  by  the  equations 

q  =  t  +  \p2-a2  (218) 

r  --  |p<  -  2a6  (219) 

s=-{f2-6z  (220) 

Elimination  of  a  and  b  leads  to  the  resolvent  cubic  equation 

t3  -  qtz  +  ( pr  -  As)t  +  4<js  -  p2s  r2  =  0  (221) 
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whence  a  and  6  can  be  calculated  either  with  the  equations 


a  =  Vjp2  +  t  -  q 
or  with  the  equations 

\pt  -  T 


b  = 


_  yt-r 


2a 


a  = 


2b 


b  =  y/\tz  -  s 


(222) 


(223) 


Computation  of  a  and  b  is  made  with  whichever  pair  of  these  equations  gives  the  least 
rounding  error.  The  choice  of  equations  is  based  upon  the  sign  of  the  expression 


yzt  -  ±qt  +  s 

If  the  coefficients  of  the  quartic  equation  satisfy  the  relationship 

2r 
P 


<?  =  iP2  + 


(224) 


(225) 


then  a  is  zero  for  one  of  the  roots  of  the  cubic  equation.  Substitution  for  q  its 
expression  in  terms  of  p  and  r  converts  the  cubic  equation  into 


0 


t3  -  (yz  +  —'if2  +  ( pr  -  4s)f  +  -  • 

V  P  /  P 

which  can  be  factored  into  the  product 

^f  -  ^  j(f2  -  ip3t  +  | pr  -  4s)  =  0 

The  roots  of  the  quadratic  factor  are 

gP2  ±  I  v'niP4  -  2 pr  +  16s 

for  which  the  square  of  b  is  given  by  the  equation 

-  s  =  iigP4  “  gPr  ±  hpZ  ^iP4  ~  2pr  +  16s 
If  the  coefficients  are  related  by  the  equation 


r 

s  =  ~5 

V 


(226) 

(227) 

(22B) 

(229) 

(230) 


then  the  radicand  is  a  perfect  square.  For  larger  values  of  s  and  a  positive  radical 
the  square  of  6  is  positive.  The  roots  of  the  quartic  equation  are  computed  with  the 
highest  real  root  of  the  cubic  equation. 

The  term  in  23  is  eliminated  from  the  quartic  equation  by  the  substitution 

z  =  <--  (231) 

4 

Then  t  is  a  solution  of  the  equation 

f4  +  af2  +  bt  +  c  -  0  (232) 


1 
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where  the  coefficients  are  defined  by  the  equations 

a  -  -  ip2  +  q  (233) 

b  =  +  sP3  “  |P9  +  r  (234) 

c  =  -  afsP*  +  ~hPB(I  -  ipr  +  s  (235) 

In  Euler's  method,  the  quartic  equation  in  t  is  factored  into  the  product 

(/  +  \  L  +  \  m  +vn)(/  +  vi  \'m  \'n)(t  -  \'l  +  Vm  ~  \/n)(t  -  Vi  -  Vm  +  Vn)  =  0  (236) 

Multiplication  of  the  factors  leads  to  the  equation 

t*  -  2(1  -t  m  +  n)f 2  +  8  v  Imn  t  +  (/  +  m  +  n)2  -  4 (mn  +  nl  +  im)  =  0  (237) 

Comparison  of  coefficients  shows  that  the  constants  i,  m,  n  are  solutions  of  the  equations 

i  +  m  +  n  =  -  |a  (238) 

mn  +  nl  +  Im  =  ^j(a2  -  4c)  (239) 

Vlmn  ---•  |6  (240) 

Thus  the  constants  l,  m,  n  are  the  roots  of  the  equation 

fc3  +  {akz  *•  /g(a2  -  4c)A:  -  ^,6*  =  0  (241) 

The  roots  of  the  quartic  equation  are  derived  directly  from  the  roots  of  the  resolvent 
cubic  equation  by  the  extraction  of  complex  square  roots. 

Programming 

SUBROUTINE  Q!CRr  (AC,  R.’) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  ROOTS  OF  QUARTIC  POLYNOMIAL 

t******************************************************************************************** 

The  real  coefficients  of  the  polynomial  are  given  in  the  array  AC  in  ascending  order. 
The  roots  of  the  polynomial  are  computed  with  the  Euler  method.  Calls  are  made  to 
Subroutine  C.BCRT.  The  real  and  imaginary  parts  of  the  roots  are  stored  in  alternate 
addresses  of  array  R.\ 


COMPLEX  ROOTS  OF  AN  ANALYTIC  FUNCTION 

Analysis 

Let  z  and  ui  be  complex  variables  which  are  expressed  in  terms  of  their  real  and 
imaginary  parts  by  the  equations 

z  x  +  iy  w  =  u  +  iv  (242) 

Let  w  be  an  analytic  function  of  c  as  expressed  by  the  equation 

w=f(z)  (243) 
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The  derivative  of  w  is  given  by  the  equation 

/  du  .  dv  \  /  dv  .  du  \ 

dw_\dx+%dx)  X  +  [dy  1  dy  )  i  dy 

-  '  ,  ,  (244) 

dx  dx  +  i  dy 

The  derivative  is  independent  of  the  direction  only  if  the  variables  satisfy  the 
Cauchy-Riemann  equations 


■  du\  w 
XJy)tdy 


du  dv 

dy  dx 


du  dv 
dx  dy 

The  modulus  of  w  is  given  by  the  equation 

[ui|  =  ( ww *)*  =  Vu2  + u* 

and  the  gradient  of  |io|  is  given  by  the  equation 

d\ui\  .  3|u>| 
V|iu[  =  — -  +  i  — — 

dx  dy 

Differentiation  and  substitution  lead  to  the  equation 

|V|ta||2  w 

VM  =  -  — r-  - 

u>  dui 


Double  differentiation  in  the  Cauchy-Riemann  equations  shows  that  u.  v  satisfy  Laplace's 
equation  as  expressed  by  the  equations 


92u  dzu 
dxz  +  9y2 

whence  the  Laplacian  of  |ia|  is  given  by  the  equation 

dz\w\  ^  dz\w\  |V|ui||2 

dxz  dy 2  |ti;| 


dzv  dzv 
dx*  dy 2 


has  the 

matrix 

d2M 

a2|ui| 

dx 2 

dxdy 

a2M 

d2|ui| 

dxdy 

dy* 

The  product  of  the  characteristic  roots  of  the  matrix  is  the  determinant 

a VI  a8M  /aVl\2  .  . 

dx 2  dy 2  \dxdy) 

At  a  point  of  zero  gradient  the  determinant  is  zero  or  negative  and  the  point  is  a 
saddle  point  unless  it  is  a  root.  The  only  minima  in  |ui|  occur  at  the  roots  of  |ui|.  The 
value  of  | in |  in  a  finite  region  is  a  maximum  only  on  the  perimeter  of  the  region  Other 
analyses  of  the  minima  and  the  maxima  are  to  be  found  in  various  texts. 


Once  a  root  has  been  established  it  must  be  eliminated  from  further  consideration. 
Let  p(z)  be  defined  by  the  equation 


Pi2)  =  I!  (2  -  <*t) 

<-1 


(253) 


where  a,,  ',at  are  the  first  k  roots.  If  10(2)  is  defined  by  the  equation 

fi2) 


w(z)  = 


p(z) 


(254) 


then  the  roots  of  w(z)  are  only  those  roots  of  f(z)  which  have  not  yet  been  computed. 
The  first  derivative  is  given  by  the  equation 


dzw  f 
dzz 


- 2- 

P  P 


dw 

=  r 

_  / 

y  1  ] 

dz 

P 

P 

(2  -  a4)_ 

is 

given 

by 

the  equation 

-  k 

y 

1 

1 

+  f 

[f  J  - 

_L_  1 

-<=1 

(2  - 

a*)J 

P 

liM  (2  -  “t). 

P 

"1  (2  -  “i)2. 

(255) 


(256) 


where  /'  and  f"  are  the  derivatives  of  /. 

During  hunting,  the  displacement  to  an  intercept  on  the  complex  plane  is  computed 
in  accordance  with  the  equation 


w 

A  2  =  -  — 
dw 


(257) 


dz 


and  the  position  2'  of  the  intercept  is  computed  in  accordance  with  the  equation 

2'  =  z  +  Az "  (258) 


For  each  displacement  Az  of  the  position  z  there  is  a  displacement  Az'  of  the  intercept  z\ 
Hunting  is  continued  as  long  as  |Az'|  is  more  than  !Az|. 

The  local  expansion  of  w  in  terms  of  z  is  given  to  second  order  by  the  equation 


dw  ,  .  ,  dzw  , 

w  +  Aui  =  w  +  -  (Az)  +  5  — j  (Az)2 

dz  dz 2 


(259) 


The  positions  of  two  nearest  roots  are  computed  when  the  equation 

w  +  Au>  =  0 


(260) 


is  solved  with  the  aid  of  the  quadratic  formula.  A  displacement  is  computed  with  the 
aid  of  the  equation 


Az 


Zw 


dw 
dz  * 


(dw\z 
\  \  dz  / 


-  2w 


d  1 u 


(261) 


dzz 


where  the  ±  sign  is  selected  to  be  whichever  sign  makes  |Az|  the  lesser.  Then  the 
displacement  is  reduced  to  a  step  of  length  <5  in  the  same  direction  as  expressed  by 


the  substitution 


***‘itoi  <262) 

The  local  representation  of  w  by  a  quadratic  approximation  eliminates  the  risk  of 
entrapment  at  a  saddle  point  in  |ia|. 

During  homing,  a  displacement  to  an  intercept  on  the  complex  plane  is  computed 
in  accordance  with  the  equation 


Az" =  ~  —— —  (263) 

dui 

dz 

and  z  is  replaced  in  accordance  with  the  substitution 

2  ■+  z'  =  z  +  Az"  (264) 

Homing  is  continued  as  long  as  \Az"\  is  less  than  |Az|.  A  tolerance  e  distinquishes 
between  disturbance  at  a  saddle  point  and  rounding  error  at  a  root.  The  position  z 
is  interpreted  as  a  saddle  point  if 

|Az"|  §  \Az\  >  e  (265) 

while  the  position  2  is  interpreted  as  a  root  if 

|Az|  S  |Az"|  <  e  (266) 

Regardless  of  the  value  of  e,  a  root  is  determined  to  the  full  accuracy  to  which  the 
functions  can  be  computed. 

When  the  point  of  computation  is  near  enough  to  a  multiple  root,  the  displacement 
in  each  cycle  is  a  constant  fraction  of  the  distance  to  go  to  the  root.  The  distance 
to  go  is  the  sum  of  an  infinite  geometric  series.  Displacement  to  the  root  in  a  single 
cycle  is  achieved  when  the  displacement  to  the  intercept  is  extrapolated  to  the  interval 


|Az|  „ 

lAz|  -  | Az  "I 


(267) 


which  expresses  the  summation  of  the  geometric  series. 

Programming 

SUBROUTINE  ENCTN  (MO,  AZ,  FT,  DF.  SF) 

t******************************************************************************************** 

EXTERNAL  SUBROUTINE  FOR  FUNCTION  AND  DERIVATIVES 
********************************************************************************************* 

The  mode  of  operation  MO  is  0  for  functions,  1  for  first  derivatives,  and  2  for  second 
derivatives.  The  argument  z  is  given  in  the  2-array  AZ.  The  function  /(z)  is  stored  in 
the  2  array  FF,  the  first  derivative  df/dz  is  stored  in  the  2-array  DF,  and  the  second 
derivative  dzf/dzz  is  stored  in  the  2- array  SF. 
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SUBROUTINE  CXRT  (CD,  CE,  AZ,  NA,  AA,  FNCTN) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  COMPLEX  ROOT  DETERMINATION 
********************************************************************************************* 

The  step  length  6  for  hunting  is  given  in  argument  CD,  and  the  tolerance  e  for  homing 
is  given  in  argument  CE.  The  initial  position  z  is  given  in  the  2-array  AZ.  The  number 
of  roots  n  is  given  in  argument  NA.  The  region  of  a  root  is  located  by  a  hunting 
procedure  with  steps  of  fixed  length,  then  the  root  is  determined  by  a  homing  procedure 
with  Newton-Raphson  iteration.  References  are  made  to  the  external  subroutine  FNCTN 
to  obtain  values  for  f(z),f'(z),f"(z ).  The  real  and  imaginary  Darts  of  the  roots  are 
stored  in  alternate  addresses  of  the  2n-array  AA. 


CONTOURS  AMONG  DATA 


The  General  Purpose  Contouring  Program23  or  GPCP  of  CalComp  works  with  a  random 
array  of  elevations.  Through  each  datum  among  the  random  data  is  passed  a  plane 
with  an  orientation  which  is  the  weighted  average  of  the  directions  toward  the  nearest 
neighbors  among  the  data.  The  weight  function  for  the  average  has  a  bell-shaped 
distribution.  The  random  array  is  converted  into  a  regular  grid  by  an  application  of 
the  weight  function  to  the  heights  of  planes  from  neighbors  nearest  to  a  grid  p  >int. 

The  elevation  in  the  interior  of  a  square  is  expressed  as  the  sum  of  products  of 
functions  which  conform  to  the  elevation  and  gradient  at  the  four  corners  of  the 
square.  For  interpolation  with  respect  to  x  the  functions  ip(x)  are  listed  for  a  unit 
square  in  the  following  table. 


v(x)  <p(0) 

(1  -  *)2(1  +  Zx)  1 

+  x2(3  -  Zx)  0 

+  x(  1  -  ar)2  0 

-*2(l-x)  0 


<p(l)  'P’(O)  <P‘(  1) 

0  0  0 

1  0  0 

0  1  0 

0  0  1 


There  is  an  analogous  table  for  interpolation  with  respect  to  y  The  products  of  the 
first  two  of  one  set  with  all  of  the  other  set  provide  enough  functions  to  meet  the 
twelve  conditions  of  elevation  and  gradient  at  the  four  corners.  Continuity  of  elevation 
and  slope  across  each  side  of  the  square  is  guaranteed  insofar  as  the  elevation  and 
slope  along  the  side  are  determined  by  the  elevations  and  gradients  at  the  ends  of 
the  side. 

The  square  is  subdivided  into  a  subgrid,  and  elevations  are  computed  at  grid  points 
with  the  cubic  approximation  for  elevation  in  the  interior.  Contours  are  traced  through 
the  subgrid  with  linear  interpolation  within  each  square  of  the  subgrid. 

The  ultimate  contouring  system  for  topographic  applications  has  been  developed 
by  Junkins  and  Jancaitis24,25.  The  Contouring  Via  the  Surface  Averaging  Concept  or 
CONSAC  system  works  with  a  rectangular  array  of  elevations  over  a  uniformly  spaced 
grid  of  rectangular  coordinates.  The  elevation  z  is  expressed  at  a  local  center  of 
approximation  by  the  equation 

z  =  c0  +  Cj  x  +  czy  +  c3xy  -  f(x,  y)  (268) 
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where  x,  y  are  Cartesian  coordinates  with  origin  at  the  center  of  approximation.  The 
coefficients  in  the  local  approximation  are  derived  by  sequential  least  squares  from 
a  local  array  of  data  which  is  centered  over  the  center  of  approximation.  Rational 
weight  factors  are  applied  to  the  data  during  the  preliminary  fitting  to  establish  the 
local  approximation  for  each  corner  of  a  square  whose  lower  left  corner  is  at  the 
origin. 

The  local  approximations  at  the  corners  of  the  square  are  blended  in  the  interior 
of  the  square.  Rational  weight  factors  are  applied  to  the  local  approximations  during 
the  final  fitting  to  establish  the  interior  approximation  for  the  whole  square.  For  a 
center  of  approximation  at  the  origin  of  coordinates  the  weight  factor  w(x ,  y)  is  given 
by  the  equation 


,  v  _  (l-»)a(l-v) 

K  ■y>  (*•+(! -*)«»»«+(!  -y)2| 


(269) 


while  the  partial  derivatives  of  the  weight  factor  are  given  by  the  equations 


dw 

2x(l  -  x) ( 1  -  y)2 

(270) 

dx 

jx2  +  (l  -x)2!2iy2  +  (l  -y)2j 

dw 

(1  -  x)22y(  1  -  y) 

(271) 

dy 

(x2  +  (1  -  x)2jjy2  +  (1  -  y)2)2 

Thus  the  value  of  w  is  unity  at  the  origin  but  is  zero  at  the  other  three  corners  of 
a  unit  square,  and  the  derivatives  of  w  are  zero  at  all  four  corners  of  the  unit  square. 
The  weight  factor  satisfies  the  equation 

w(x.  y)  +  w(  1  •-  x,  y)  ui(x,  1  -  y)  ±  w(  1  -  x,  1  y)  =  1  (272) 

throughout  the  interior  of  the  square. 

The  elevation  in  the  interior  of  the  square  is  given  by  the  equation 

z  -  wji  +  wz f2  +  Wsfj  +  uiJt  =  F(x,  y)  (273) 

where  wt,  wz,  w3.  wt  are  the  weight  factors  for  each  of  the  four  corners  and /|,/g,/3,/4 
are  the  local  approximations  for  each  of  the  four  corners.  Since  the  weight  factors 
have  quadratic  numerators  and  denominators  while  the  local  approximations  are 
linear,  the  final  approximation  for  the  interior  of  the  square  is  the  ratio  between  a 
cubic  polynomial  and  a  quadratic  polynomial. 

Differentiation  with  respect  to  x  along  a  line  of  constant  y  and  reduction  to  a 
common  denominator  replaces  the  numerator  of  F{x,y)  by  a  polynomial  which  is 
quartic  in  x,  while  differentiation  with  respect  to  y  along  a  line  of  constant  x  and 
reduction  to  a  common  denominator  replaces  the  numerator  of  F(x,  y)  by  a  polynomial 
which  is  quartic  in  y.  Solution  of  these  quartic  equations  gives  the  points  where  F(x,  y) 
has  minima  or  maxima  along  a  side  or  along  a  median  of  the  square.  From  the 
elevations  of  the  minima  or  the  maxima  are  derived  the  range  of  contour  levels  which 
intercept  the  sides  or  the  medians  of  the  square. 

On  a  contour  of  level  h  the  elevation  is  expressed  by  the  equation 

F(z.  y)  =  h  (274) 


Clearance  of  the  denominator  in  this  equation  gives  a  cubic  equation  for  the  contour 
line.  The  cubic  equation  is  solved  for  y  with  x  set  to  a  sequence  of  values  or  for  x 
with  y  set  to  a  sequence  of  values  to  obtain  a  sequence  of  points  on  the  contour. 
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CONTOURS  OF  FUNCTIONS 

Analysis 

Let  x,  y  be  Cartesian  coordinates  and  let  z  be  a  constant  contour  level.  Let  i,j  be 
unit  vectors  in  the  directions  of  increasing  x,  y.  Let  f(x,  y)  be  a  function  of  x,  y.  Then 
the  equation 


z^f(x.y)  (275) 

defines  implicitly  a  relation  between  x  and  y.  A  position  vector  r  is  defined  by  the 
equation 

r  -  xi  +  yj  (276) 

The  differential  change  df  in  the  function  /(x,  y)  for  the  differential  change  dr  in  the 
position  vector  r  is  given  by  the  equation 


3/  df 
df  =  —  dx  -  --dy 
J  dx  dy  * 


V/  dr 


(277) 


and  the  gradient  of  /(x,  y)  is  defined  by  the  equation 

V/=i|^+j“  (278) 

dx  dy 

For  each  trial  point  there  is  an  intercept  in  the  x,  y  plane  where  the  gradient  is 
estimated  to  be  zero.  The  displacement  Ax*,  Ay*  from  the  trial  point  to  the  point  of 
zero  gradient  is  a  solution  of  the  equations 


a2/  a2/  df 

3x2  dxdy  dx 


a2/  a2/ 

J  Ax*  +  -  2  Ay* 


dxdy  3y2 

The  solution  is  expressed  by  the  equations 

Of  d^f 

dx  3y2 


Of 

dy 


Ax*  = 


df  a2/ 

dy  dxdy 


(279) 

(280) 

(281) 


dy  8x2 

Ay  a2/  a2/ 


(282) 


a2/  a2/  /  a^M2 

3x2  3y2  \dxdy) 

df  d 7  df  a2/ 

3x  dxdy 

(diJ\Z 

dxz  dy2  \dxdy) 

The  coordinates  x*,y*  of  the  point  of  zero  gradient  are  given  by  the  equations 

x*  -  x  +  Ax*  y*  -=  y  +  Ay*  (283) 

The  value  f*  of  the  function  at  the  coordinates  x*,y*  is  given  by  the  equation 

r  -  f  +  r_  Ax*  +  r::  Ay*  +  ‘  ■  ~  (Ax*)2  +  Ax* Ay*  *  \  0  (Ay*)2  (284) 
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The  matrix  of  the  equations  which  determine  the  point  of  zero  gradient  is  the  matrix 
of  the  gradient  of  the  gradient  of  /.  It  is  symmetric,  and  it  has  real  characteristic 
roots  with  orthogonal  characteristic  vectors.  Its  trace  is  the  sum  of  its  characteristic 
roots,  and  its  determinant  is  the  product  of  its  characteristic  roots. 

Each  of  the  characteristic  roots  y,  v  of  the  matrix  is  given  by  the  equation 

i  /  a2e  a2e\  i  f/~a2/  a2r\2  /  a*/  \2 


i/av  .  «V\  i 

2\3x2  ~  3y2/  2 


p/. 

d2/V.4/ 

<M 

> 

I\dx2 

dyz)  \ 

,  dxdy ! 

where  the  choice  of  root  is  determined  by  the  ±  sign.  The  characteristic  unit  vectors 
m,  n  are  given  by  the  equation 


!  ,  Va2/ 

a2/\ 

f  i  J 

'a2/ 

a2/\2 

J  \2  /'2 

[  *  2 \ dx2 

aW 

2  \'\ 

Kdxz 

ay2/ 

V dxdy j  1 

(/a2/ 

a2/' 

,  2 

4.4( 

'  a2/  \2) 

1 

4 

(lax2 

ay2, 

1  4\ 

1 3x3y  /  ) 

UI  (dZl 

-a2/u 

ijl 

-a2/ 

dZf\Z  , 

4/  av )  n 

(  2\dx2 

ay2/ 

2  Vi 

,  3x2 

dyz ) 

l dxdy /  ) 

{(dzf 

a2/' 

'  \z) 

1 

(\3x2 

ay2, 

)  t4( 

,dxdy)  ) 

where  the  choice  of  vector  is  determined  by  the  ±  sign.  The  gradient  of  the  gradient 
of  f  is  given  by  the  equation 

VV/  -  /umm  i  imn  (287) 

If  the  roots  have  opposite  signs,  the  point  of  zero  gradient  is  a  saddle  point,  and  the 
plane  of  z  intersects  the  surface  of  /.  If  the  roots  have  the  same  signs,  the  point  of 
zero  gradient  is  a  minimum  or  a  maximum.  If  the  signs  of  the  roots  then  are  the 
same  as  the  sign  of  /  z,  the  plane  of  z  docs  not  intersect  the  surface  of  f 

Let  p  be  the  position  of  the  trial  point  relative  to  the  point  of  zero  gradient  as 
expressed  by  the  equation 


Let  u.v  be  Cartesian  coordinates  with  unit  vectors  m.n  in  the  directions  of  increasing 
u,  u.  The  position  vector  p  of  the  trial  point  is  given  by  the  equation 

p  --  um  vn  (289) 

The  coordinates  u,  v  are  expressed  in  terms  of  the  coordinates  x,  y  by  the  equations 

u  -  ( x  x*)i  m  t  (y  -  y*)j  m  (290) 

v  =  (x  x*)  i  n  +  (y  y*)  j  n  (291) 

Let  be  the  value  of  the  function  relative  to  the  point  of  zero  gradient  as  expressed 
by  the  equation 

V-f  r  (292) 

The  value  of  the  function  is  given  in  terms  of  the  coordinates  by  the  equation 

J  l2vvz  (293) 

If  the  characteristic  roots  have  opposite  signs  the  contours  of  constant  <p  are  hyperbolas 
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and  if  the  roots  both  have  the  same  sign  as  p  the  contours  of  constant  p  are  ellipses. 

The  trace  of  the  matrix  of  the  equations  which  determine  the  point  of  zero  gradient 
is  given  by  the  expression 


dx 2 


8V 

3y 2 


and  the  determinant  of  the  matrix  is  given  by  the  expression 

aV  a2/  /  a2/  \2 


dx*  dyz 


\ dxdy  ) 


(294) 


(295) 


If  the  determinant  is  zero  there  is  no  single  point  of  zero  gradient.  The  partial 
derivatives  are  related  in  accordance  with  the  equation 


a2/ 

dxdy 


i  i 

\2/ 


*  lax2 


n 


dy'‘ 


and  the  function  is  given  for  any  Ax*.  Ay*  by  the  equation 


(296) 


/*-=/-  2 


a  a//av\2 
*  ay  Uy2j 

a2 

dx 

f  /32P 
r lay2/ 

/  _  a2/ 

2  ay2 

%  V/aV\s 

ay  lax2/ 

a2 

ax 

r  /  a  2/  \ 

r  lax2, 

/  a2/ 

2  ay2 

3y  \3y  / 

a*/\*  ,/a7\* 


3y‘ 


r:»  Ax**  ITS  Ay* 


ax2 


a2/ 

3y2 


(0) 


Ax*  ± 


\dx‘ 


a2/\s  . 

— t)  Ay* 


ay 


(297) 


The  contours  of  constant  function  are  parabolas.  The  nearest  point  on  the  axis  of  the 
parabolas  is  reached  when  Ax*,  Ay*  are  given  by  the  equations 


Ax*  =  - 


a/  /  av\*  of  (dy_\z 

dx  [dx1/  dy  \3y2/ 


Ay*  =  * 


V/a*/ 

3x  \ 3x2 


/aV  ,  aVv 

lax2  ay2/ 

v/»v\* 

/  dy\dyzj 


(3) 


lax2 


av 

ay2 


m 


(298) 


(299) 


Both  terms  of  the  trace  have  the  same  sign  and  their  square  roots  are  positive  or 
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imaginary.  Let  X,  v  be  defined  by  the  equations 


X  = 


a/  /a2/y>  ^  df  (dy  y 

_ dy  lax2/ 

(dl ,  «v\* 

l~3x2  dy 2/ 


a2/  a2/ 

Si2  dy2 

and  let  m,  n  be  defined  by  the  equations 

i 


(300) 


(301) 


n  ± 


/a2  A* 

(*fy 

lay2/ _ 

lax2/ 

'  7  d*f\ilT 

/  a2/ 

<x2  ~  ay2/ 

lax2  3y2, 

(a*fy 

/a2/ J 

\dxz  I 

,  >  ■*- 

lay2/ 

/a2/  4 

a2/\s 

/a8/ 

a2/\ 

laz2 

ay2/ 

l  az2 

ay2/ 

(302) 


(303) 


Let  the  coordinates  u,v  be  derived  from  the  coordinates  x.y  in  accordance  with  the 
equations 


u  ( x  x*)i  m  -•  (y  -  y’)j  m 
V  (z  x*)i  n  4  (y  -■  y*)j  n 
Then  the  function  p  is  given  by  the  equation 

Ifi  Xu  4  lzUVZ 


(304) 

(305) 

(306) 


The  contours  of  constant  ip  are  parabolas  with  axis  in  the  direction  of  m. 

During  hunting  perpendicular  to  a  contour  a  displacement  to  an  intercept  on  the 
x.  y  plane  is  computed  in  accordance  with  the  equation 


Ar"  a  :_V,  V/  (307) 

and  the  position  of  the  intercept  is  computed  in  accordance  with  the  equations 


-l—  Vf 
;v/2  1 


X  -  Ax" 


y'  =  y  +  Ay” 


(306) 


The  increments  Ax,  Ay  are  replaced  by  a  step  of  length  6  in  the  direction  of  Ar"  as 
expressed  by  the  equation 


Ar  -  6 


Ar" 

!Ar"' 


(309) 


unless  a  point  of  zero  gradient  is  detected  within  that  circle  which  is  defined  by  the 
equation 


(z*  z)2  ‘  (y*  y)2 


(310) 
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The  characteristic  roots  and  vectors  of  the  matrix  of  Wf  are  analysed  to  determine 
whether  the  plane  at  z  can  intersect  the  surface  of  /.  Computation  is  terminated  if 
there  can  be  no  intersection.  Otherwise  the  trial  point  is  displaced  toward  the  nearest 
point  of  intersection.  Hunting  is  continued  as  long  as  |Ar'|  is  more  than  |Ar|,  then 
homing  is  started  when  |Ar’|  becomes  less  than  I Arj. 

During  hunting  parallel  to  a  contour  the  steps  of  the  trial  point  are  expressed  by 
quadratic  curves  which  are  osculatory  to  the  contour. 

If  /u,  v  have  opposite  signs,  then  the  contour  of  constant  <p  is  an  hyperbola.  Let  the 
sign  of  be  the  same  as  the  sign  of  <p  so  that  u,  v  can  be  expressed  in  terms  of  a 
parameter  17  by  the  equations 
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u  =  ^  — y  cosh  77 

(311) 

V  2  sinh  T] 

1  i*.l 

(312) 

Increments  Au,  Au  in  the  coordinates  are  given  in  terms  of  the  increment  At]  in 
parameter  by  the  equations 


tZ<pXz 

A u  -  j  sinh  77  Atj 

(313) 

!  2y?  1 1 

Au  -  1  cosh  77  A77 

1  v  | 

(314) 

The  hyperbolic  functions  may  be  expanded  with  the  aid  of  the  addition 
may  be  expressed  in  terms  of  the  coordinates  to  give  the  equations 

formula  and 

!  u  |  2 

u  +  Au  ■-  u  cosh(A77)  +  —  |  u  sinh(A7j) 

Mi 

(315) 

\  Li\  2 

v  +  Au  -  u  cosh(A77)  +  -  u  sinh(A77) 

1  v  1 

(316) 

which  are  valid  for  arbitrary  increments,  and  the  equations 

1 

■  V  I  2 

Au  -  1  1  v  At] 

1  M  i 

(317) 

/U  !  2 

Au  -•  -  i  u  At] 

(318) 

I  v  i 


which  arc  valid  for  small  increments  The  increments  Au,  Av  are  initially  in  the  direction 
of  the  tangent,  which  can  be  derived  directly  from  the  gradient  V/.  The  stability  of 
the  computation  is  controlled  when  the  increments  are  subject  to  the  limitation 

\  (Au)2  •  (An)2  '  <5  (319) 

where  (5  is  the  maximum  length  of  step  Then  the  increment  At]  is  determined  with 


i 


17 


1 

f 


sufficient  accuracy  by  the  equation 


Ar)  =  - 


Ati 


v(A u)2  +  (Au)2  I 


u 


Au 


\  (Au)2  -  (Au)s 


(320) 


unless  this  would  make  , Ar7  more  than  in  which  case  At;  is  reduced  to  ± 

If  is  zero,  the  contour  of  constant  <p  is  a  parabola  Let  u,  u  be  expressed  in  terms 
of  a  parameter  a  by  the  equations 


u  = 


X 


X 


v  -  -a 
v 

Increments  Au,  At.  in  the  coordinates  are  given 
parameter  by  the  equations 


(321) 

(322) 

in  terms  of  the  increment  A  a  in 


Au - ct  Act 

v 

X 

Au  4  -  Act 
v 

The  coordinates  are  given  by  the  equations 


(323) 

(324) 


u  +■  Au  -  u  u  A  ct  - 


1  X, 

5  (5ct)2 

2  v 


v  +  Au  =  u 


X 

-  Act 
v 


(325) 

(326) 


which  are  valid  for  arbitrary  increments  and  the  increments  are  given  by  the  equations 


Au  -  •-  v  Act  (327) 

X 

Au  -  +  -  Act  (328) 

v 


which  are  valid  for  small  increments.  The  increments  Au,  Au  are  initially  in  the  direction 
of  the  tangent,  and  the  length  of  step  is  limited  to  <5.  Then  the  increment  Act  is 
determined  with  sufficient  accuracy  by  the  equation 


Act 


_ <5 _ 

Au  X  Au 

V(Au)2  +  (Au)2  v  V(Au)2  +  (Au)2 


(329) 
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unless  this  would  make  'Act|  more  than  in  which  case  Act  is  reduced  to  ± 

If  n,  v  both  have  the  same  sign  as  <p,  then  the  contour  of  constant  <p  is  an  ellipse. 
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unless  the  initial  point  is  detected  ahead  on  the  path  of  hunting.  When  the  position 
x0,  y0  of  initiation  lies  within  that  zone  which  is  defined  by  the  criteria 


(x0  x)2  +  (y0  y)2  <  (Ax)2  +  (At/)2  (342) 

l(*o  -  x)Ay  -  (y0  y)Ax|  <  c  V (Ax)2  +  (Ay)2  (343) 

where  e  is  a  specified  tolerance,  then  the  coordinates  x,  y  are  replaced  in  accordance 
with  the  substitutions 

*  *  *o  y  ’Vo  (344) 

and  the  computation  is  terminated.  Otherwise  the  hunting  procedure  is  replaced  by 
the  homing  procedure. 

During  homing  the  coordinates  x.  y  are  replaced  in  accordance  with  the  substitutions 
x  >  x'  y  >  y'  (345) 

until  Ax,  Ay  satisfy  the  criterion 


Ar  ■  c  (346) 

and  the  displacements  of  the  trial  point  have  been  reduced  to  the  level  of  rounding 
error. 

If  the  trial  point  is  stepped  out  of  bounds,  it  is  reset  on  the  boundary  If  x  is 
outside  the  boundary  it  is  reset  onto  the  boundary  and  y  is  adjusted  in  linear  proportion, 
or  if  y  is  outside  the  boundary  it  is  reset  onto  the  boundary  and  x  is  adjusted  in 
linear  proportion.  The  computation  is  terminated  if  the  trial  point  is  trapped  in  a 
corner.  If  x  already  is  on  a  boundary,  then  partial  derivatives  with  respect  to  x  are 
set  equal  to  zero,  or  if  y  already  is  on  a  boundary,  then  partial  derivatives  with  respect 
to  y  are  set  equal  to  zero.  Intersection  of  the  contour  with  the  boundary  is  established 
by  iteration 

After  Newton  Raphson  iteration  the  function  and  its  derivatives  are  known  at  each 
trial  point.  Two  consecutive  trial  points  are  centers  of  expansion  from  each  of  which 
the  contours  arc  extended  halfway  to  the  other.  Insofar  as  the  extensions  arc  accurate 
only  to  second  order,  they  do  not  quite  meet  in  the  middie.  The  extensions  are  adjusted 
in  proportion  to  the  cube  of  the  distance  from  their  centers  of  expansion  in  order  to 
bring  them  into  coincidence  midway  between  the  trial  points 

Programming 

SUBROUTINE  fNCHJ  (MO,  AX,  A>.  ft,  tv,  \>.  «.  • ,  ■'■) 

EXTERNAL  SUBROUTINE  FOR  SURFACE  ELEVATION  AND  DERIVATIVES 

The  mode  of  operation  VO  is  0  for  functions,  I  for  first  derivatives,  and  2  for  second 
derivatives.  The  coordinates  x,  y  are  given  in  the  arguments  •'•-.a-.  The  elevation  /  is 
stored  in  the  function  :  I  .  The  first  derivatives  df  dx.  df  dy  arc  stored  in  the  functions 
.  the  second  derivatives  d2/  dx2.  d2/  Sidy.  d2/  dy2  are  stored  in  the  functions 
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SUBROUTINE  CNTCRV  (CD.  CF,  CF.  CA,  AX.  AY.  A/,  MO,  ENCTN) 

**  +  +  **  +  *********  +  +  +  *  +  ****  +  *  +  +  +  +  +  +  ****  +  *  +  +  +  *if  +  ***-irm-t.**  +  +  *m,m**  +  +  *  +  ii****tm*  +  mm  +  mtm 

FORTRAN  SUBROUTINE  FOR  CONTOUR  CURVE  CONSTRUCTION 

***«*******************************************4********************************  ************* 

The  step  length  6  for  hunting  is  given  in  argument  CD,  and  the  tolerance  e  for  homing 
is  given  in  argument  CF.  The  edges  of  the  plotter  field  are  given  in  the  4- array  CF  in 
the  order  left,  right,  upper,  lower,  and  the  edges  of  the  contoured  area  are  given  in 
the  4  array  CA  in  the  order  left,  right,  upper,  lower.  The  initial  coordinates  x„  and  Vo 
are  given  in  the  arguments  AX  and  AY.  The  elevation  of  the  contour  is  given  in  the 
argument  A.'.  Datum  points  on  the  contour  are  written  on  tape  file  NO.  The  elevation 
of  the  surface  is  computed  by  reference  to  the  external  subroutine  FNCTN. 


DISCUSSION 

Experiments  have  been  performed  on  the  computer  to  determine  the  levels  of 
rounding  error  in  power  polynomial  approximation.  Values  of  polynomials  at  nodes 
and  at  antinodes  and  coefficients  of  polynomials  have  been  compared  for  a  few  orders. 
In  all  of  the  experiments  the  abscissae  ranged  from  1  to  fl. 

For  equal  spacing  the  values  of  the  Lagrange  polynomials  are  greatest  for  those 
polynomials  which  are  unity  at  abscissae  near  the  middle  of  the  range  of  abscissae. 
The  values  at  antinodes  are  larger  than  unity  but  less  than  the  maximum  values  at 
endpoints. 

For  Chcbyshcv  spacing  the  values  of  the  Lagrange  polynomials  are  greatest  for 
those  polynomials  which  arc  unity  at  abscissae  near  the  ends  of  the  range  of  abscissae. 
The  values  at  antinodes  arc  everywhere  less  than  unity  for  all  polynomials  but  the 
values  at  endpoints  are  greater  than  unity. 

The  increase  in  endpoint  values  with  order  is  illustrated  by  the  following  table. 


Endpoint  Maxima 


Number 

Equal 

Chcbyshcv 

Degree 

of  Data 

Spacing 

Spacing 

1  1 

8.6-  10' 

1.27 

10 

17 

3.5-  103 

1.27 

16 

21 

4.5'  10* 

1.27 

20 

25 

6.1  •  105 

1.27 

24 

Thus  any  error  in  an  ordinate  near  the  middle  of  the  range  of  abscissae  is  amplified 
into  wild  wiggles  near  the  ends  of  the  range  when  the  abscissae  arc  equally  spaced 
The  efTcet  of  rounding  error  is  related  to  the  disparities  between  the  magnitudes 
of  terms  and  the  values  of  polynomials.  The  increase  in  the  coefficients  of  terms  of 


Lagrange  polynomials  with  order  is  illustrated  by  the  following  table. 

Maximum  Coefficients 


Number 

Equal 

Chebyshev 

Degree 

of  Data 

Spacing 

Spacing 

1 1 

3.2x  103 

2.6x  102 

10 

17 

1.4x10® 

2.9x10* 

16 

21 

1.0/10® 

7.4x10® 

20 

25 

6.8x10® 

1.9/107 

24 

The  increase  in  the  coefficients  of  terms  of  orthonormal  polynomials 
illustrated  by  the  following  table. 

Maximum  Coefficients 

with  order 

Number 

Equal 

Chebyshev 

Degree 

of  Data 

Spacing 

Spacing 

11 

5.3  <103 

5.5xl02 

10 

17 

2.7x10® 

7.3x10* 

16 

21 

2.2xl07 

2.0x10® 

20 

25 

1.4xl010 

5.6xl07 

24 

The  rounding  error  for  abscissae  near  the  ends  of  the  range  of  abscissae  is  on  the 
order  of  the  product  of  the  maximum  coefficient  and  Z'N  for  A'- bit  arithmetic.  Thus 
expansion  into  coefficients  becomes  rapidly  more  vulnerable  to  rounding  error  with 
increasing  order. 

If  the  abscissae  are  proportional  to  the  cosines  of  equally  spaced  angles,  then  the 
trigonometric  functions  of  multiples  of  the  angles  are  orthogonal  with  respect  to 
summation,  and  the  orthonormal  polynomials  are  proportional  to  the  trigonometric 
functions.  The  three  term  recurrence  for  the  polynomials  is  just  the  addition  theorem 
for  the  trigonometric  functions.  Experiments  on  the  CDC  6600  with  exact  recurrence 
separate  the  effects  of  rounding  error  in  the  three  term  recurrence  from  the  effects 
of  rounding  error  in  the  evaluation  of  coefficients.  That  the  rounding  error  in  the 
recurrence  increases  only  very  slowly  with  increasing  order  is  indicated  by  the  following 
table. 


Number  of  Data 

Maximum  Error 

Degree 

11 

1.3xl0“13 

10 

17 

2. Ox  10"13 

16 

21 

3.2/  10'13 

20 

25 

6.9x  10-13 

24 

where  the  maximum  error  is  estimated  from  the  residual  at  the  nth  abscissa  after 
the  (n  t  l)th  cycle  of  recurrence. 

There  seem  to  be  as  many  subroutines  for  the  EFT  as  there  are  users.  Only  three 
of  the  existing  subroutines  will  be  considered  in  the  present  discussion.  A  subroutine 
NLOGN  has  been  published  by  Robinson11.  It  is  used  extensively.  It  is  short  and  moderately 
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fast.  A  subroutine  FFT  has  been  published  by  Brigham12.  It  is  relatively  straightforward. 
It  is  short,  but  quite  slow.  An  improved  version  CXFFT  has  been  prepared  for  the  present 
program.  One  major  subroutine  is  FFTP,  which  is  distributed  by  1MSL.  It  is  very  long 
but  very  fast.  It  is  not  limited  to  orders  which  are  a  power  of  3.  It  factors  any  order 
into  prime  factors  and  applies  the  original  formulation  of  Cooley  and  Tukey.  The 
ultimate  subroutine  is  HARM,  which  is  distributed  by  IBM.  It  is  very  long  but  very  fast. 
It  is  three  dimensional.  A  set  of  three  subroutines  has  been  prepared  for  the  present 
program.  They  can  be  used  in  computing  loops  to  perform  the  same  functions  as 
HARM.  The  computing  loops  are  only  half  as  fast,  however,  because  they  do  not  take 
every  possible  shortcut  in  computation.  The  merit  of  each  subroutine  depends  upon 
the  length  of  program,  the  storage  requirement  for  data,  and  the  speed  of  computation. 
Any  increase  in  speed  of  computation  is  achieved  at  the  expense  of  an  increase  in 
storage  requirement. 

In  the  IMSL,  subroutine  ZANLYT  applies  the  Muller  method,  and  subroutine  i'CPOLY 
applies  the  method  of  Jenkins  and  Traub.  In  the  present  system,  CXRT  follows  the 
hunting  and  homing  method.  The  length  of  ZANLYT  is  a  few  percent  less  than  the  length 
of  CXRT,  but  the  length  of  ZCPOLY  is  more  than  three  times  the  lengths  of  the  other 
two  subroutines.  Comparative  tests  with  tenth-degree  polynomials  showed  that  ZCDOL« 
takes  a  few  percent  less  time  than  the  other  two  subroutines  For  roots  equally  spaced 
along  the  real  axis,  all  three  subroutines  gave  full  machine  accuracy,  but  for  roots 
equally  spaced  around  a  circle  of  unit  radius,  ZA\;.»"  and  ,'CPO.  ■  gave  two  to  three 
digits  less  accuracy  than  CXRT.  All  three  subroutines  gave  double  roots  with  only  half 
the  accuracy  of  single  roots. 


CONCLUSION 

Although  Lagrange  and  orthonormal  polynomials  of  high  degree  should  not  be 
expanded  in  coefficient  form,  they  can  be  evaluated  accurately  to  any  order  by 
recurrence  relations.  Insofar  as  computations  of  roots  and  contours  by  hunting  and 
homing  procedures  are  applied  to  unreduced  functions,  they  give  high  accuracy 
determinations. 


I 
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